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1  INTRODUCTION 

POINT  is  the  acronym  for  a computer  program  for  orbit  integration.  The 
program,  which  numerically  integrates  the  equations  of  motion  of  an  earth  satel- 
lite in  a highly  eccentric  orbit,  was  originally  conceived  as  a component  of  the 
computer  software  required  for  the  analysis  of  the  launch  phase  of  a synchronous 
satellite  mission.  The  program  was  written  with  a view  to  more  general  applica- 
tion, however,  and  will  be  useful  in  pre-launch  analyses  of  other  satellite 
projects.  A shortened  version  has  already  been  used  in  this  role  to  assess  the 
orbital  stability  of  the  proposed  Astronomical  Roentgen  Observatory  (ASRO) * . 

Given  a set  of  initial  conditions,  the  program  generates  ephemerides  of  an 
earth-orbiting  satellite.  At  the  epoch,  the  orbit  can  be  specified  either  by 
the  cartesian  components  of  the  geocentric  position  and  velocity  vectors,  by  a 
standard  set  of  launch  vehicle  injection  conditions  or  by  one  of  three  types  of 
orbital  elements.  The  ephemerides,  of  which  any  combination  may  be  obtained, 
consist  of  the  cartesian  components  of  geocentric  position  and  velocity,  osculat- 
ing orbital  elements,  look  angles  for  a maximum  of  four  ground  stations  and  a 
table  of  apses. 

Orbit  development  is  obtained  by  the  numerical  integration  of  the  equations 
of  motion  as  formulated  in  the  Cowell  method.  Large  variations  occur  in  the 
satellite’s  velocity  and  in  the  perturbing  forces  around  a highly  eccentric 
orbit.  At  perigee  where  the  velocity  is  high,  and  the  forces  are  changing 
rapidly,  a small  integration  step  length  is  required.  Use  of  this  step  length 

around  the  entire  orbit  is  unnecessary  and  inefficient,  however.  Analytical 

2  ... 
step  regulation  is  therefore  used  to  obtain  the  required  accuracy  without  loss 

of  efficiency. 

Luni-solar  perturbations  caused  by  the  differential  attractions  of  the  sun 

and  moon  are  computed  as  required  during  integration.  The  sun-moon  coordinates 

are  stored  at  daily  intervals  in  a disc  file.  The  coordinates  originate  from 
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the  JPL  ephemeris  tapes  ' and  have  been  transformed  into  PROP  axes  (see 

section  4).  Although  solar  radiation  pressure  is  not  incorporated  at  present, 

the  program  is  written  in  a manner  which  will  enable  it  to  be  inserted  at  a 

future  date.  The  asphericity  of  the  earth  is  taken  into  account  approximately 

by  including  the  effects  of  zonal  harmonics  up  to  and  tesseral  harmonics  up 

to  J,  , . Atmospheric  drag  is  included  with  the  density  evaluated  from  a 
. . 4 5 

modified  form  of  the  Jacchia  1965  model  atmosphere. 


POINT  is  written  in  1900  FORTRAN  for  use  on  ICL  1900  series  computers. 

The  program  requires  approximately  20K  words  of  core  store,  and  consists  of  a 
main  program  and  33  subprograms.  The  program  units  are  listed  in  Appendix  A, 
while  the  calling  structure  is  illustrated  in  Fig.l.  A flow  chart  for  the  main 
program  is  given  in  Fig. 2.  The  subprogram  specifications  are  given  in 
Appendix  B. 

2 PROGRAM  FUNCTION 

The  program  may  take  as  input,  at  a given  epoch,  either  the  geocentric 
cartesian  components  of  the  position  and  velocity  vectors  (x,y,z,x,y,z) , or  a 
standard  set  of  launch  vehicle  injection  conditions  or  one  of  three  types  of 
orbital  elements.  The  launch  vehicle  injection  conditions  consist  of  speed, 
climb  angle,  azimuth  of  the  velocity  vector,  radial  distance,  geocentric 
latitude  and  longitude.  The  three  types  of  orbital  elements  are:  firstly,  a set 
of  osculating  elements  defined  by  the  instantaneous  position  and  velocity  vec- 
tors, and  consisting  of  semi  major  axis  a , eccentricity  e , inclination  i , 
right  ascension  of  the  ascending  node  ft  , argument  of  perigee  a>  and  mean 
anomaly  M ; secondly,  a set  of  two-line  elements  in  the  form  issued  by  the  USAF 
Space  Defense  Center  (SDC) ; and  thirdly,  a set  of  RAE  five-card  elements  pro- 
duced as  output  by  the  computer  program  PROP  . The  latter  are  mean  elements  in 
the  sense  defined  by  Kozai^,  i.e.  time  averages  per  revolution  of  the  osculating 
elements. 

The  program  generates  up  to  three  types  of  ephemerides  each  of  which  is  a 
table  of  'satellite  positions'  at  uniform  intervals  of  time.  For  the  purpose  of 
this  Report  'satellite  position'  means  either  the  cartesian  components  of  the 
geocentric  position  and  velocity  vectors,  a set  of  osculating  orbital  elements, 
or  a set  of  look  angles  for  a given  ground  station.  Look  angles  consist  of 
azimuth  A , measured  clockwise  from  true  north,  elevation  E , range  and  range 
rate.  The  program  can  also  determine  the  satellite's  position  and  time  at  an 
apse;  this  information  is  then  output  in  tabular  form.  Although  only  three  type 
of  ephemerides  can  be  printed  at  present,  it  would  be  straightforward  to 
increase  the  program's  flexibility  by  including  other  types. 

For  cartesian  components  it  is  possible  to  obtain  the  covariance  matrix  as 
an  additional  output  provided  that  the  covariance  matrix  of  the  initial  condi- 
tions is  available  as  part  of  the  input.  A word  of  caution  is  needed  here, 
however.  Since  the  propagation  of  errors  by  means  of  a covariance  matrix  is 
basically  a linear  concept,  difficulties  can  be  encountered  if  the  initial 
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errors  are  large,  e.g.  as  may  result  from  the  errors  associated  with  a set  of 
launch  vehicle  injection  conditions. 
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3 EQUATIONS  OF  MOTION  AND  INTEGRATION  PROCEDURE 

For  an  earth  satellite  the  equations  of  motion  may  be  expressed  in  vector 
form  as 

£ = ~wr1 * 3£+F  , (1) 

where  r'  is  the  position  vector  relative  to  the  earth's  centre  of  mass, 

U (=  GM)  is  the  gravitational  constant  of  the  earth,  and  F is  the  perturbing 
acceleration.  Orbit  development  is  obtained  by  the  numerical  integration  of  the 
equations  of  motion  as  formulated  in  the  Cowell  method.  Because  of  the  large 
variations  in  the  velocity  and  perturbing  accelerations  around  a highly  eccentric 
orbit,  a constant  time  interval  cannot  be  used  if  the  computation  is  to  be 
efficient.  To  overcome  this  difficulty,  the  equations  are  transformed  in  such  a 
way  that  a constant  step-length  can  be  used,  i.e.  analytical  step  regulation  is 
introduced.  Time  t is  replaced  by  the  independent  variable  s , defined  by 

dt/ds  = rk/ x , (2) 

. . 2 

where  * is  a constant.  Merson  has  discussed  the  selection  of  k and  < to 

1 / 2 

give  the  best  results  and  suggested  the  use  of  k = 3/2  and  < = u ; these 
values  are  used  in  the  program.  On  changing  to  the  independent  variable  s , 
where 

t’  = dt/ds  = p"1/2r3/2  , (3) 

the  equation  (1)  becomes 


The  equation  for  t'  can  be  differentiated  to  give 

t”  = | ( £ . r'myr)*  , (5) 

so  that  we  have  four  second-order  differential  equations  (for  x" ,y" ,z" , t") . The 
integration  is  based  on  an  eighth-order  Gauss-Jackson  second-sum  process. 
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A sixth-order  Butcher  process  is  used  to  set  up  the  difference  table  required 
before  the  second-sum  procedure  can  be  started. 

4 TIME  AND  COORDINATE  SYSTEMS 

Calendar  dates  are  reckoned  in  Modified  Julian  Days  (MJD) , which  are 
related  to  (ordinary)  Julian  Days  by  the  formula 

MJD  = JD  - 2400000.5  . 

Q 

The  coordinate  system  used  internally  is  the  one  suggested  by  Kozai  . The 
origin  0 is  at  the  earth's  centre  of  mass  and  the  Oz  axis  points  towards  the 
north  pole.  Ox  lies  in  the  plane  of  the  true  equator  of  date,  but  instead  of 
pointing  towards  the  true  equinox  of  date  it  is  directed  towards  a projection  of 
the  mean  equinox  of  the  epoch  1950.0  (MJD  33281.9234);  Oy  completes  the  right- 
handed  system  Oxyz  . 

5 UNITS  AND  CONSTANTS 
5.  1 General 

Distances  are  measured  in  kilometres  (for  input  of  a covariance  matrix, 
however,  distances  may  alternatively  be  given  in  feet).  Angles  are  held  as 
radians  inside  the  computer,  but  are  measured  in  degrees  for  both  input  and 
output.  Time  is  measured  in  seconds,  except  when  identifying  specific  instants 
or  time  intervals. 

Constants  used  in  the  program  are  assigned  only  once,  either  in  the  main 
program  or  in  the  BLOCK  DATA,  subprogram. 

5.2  Geophysical  constants 

Values  of  the  geophysical  constants  are  assigned  in  the  BLOCK  DATA  sub- 
program; current  values  are  listed  below,  E followed  by  an  integer  designating 
an  exponent.  Values  for  the  even  zonal  harmonics  have  been  taken  from  Ref. 9, 
while  those  for  the  odd  zonal  harmonics  have  been  taken  from  Ref. 10.  Values  for 
the  tesseral  harmonics  have  been  taken  from  a set  derived  by  Wagner''. 
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Ccophysical  constants 


Cons  tant 


FORTRAN 

variable 


Value 


3 - 2 

Earth's  gravitational  constant  u (km  s ) 

EMU 

398601.3 

Mean  equatorial  radius  of  earth  R (km) 

ERAD 

6378.163 

Earth's  zonal  harmonics  J 2 

EJ2 

1 082 . 637E-6 

J3 

EJ3 

-2.531E-6 

EJ4 

-1 . 6 1 9E-6 

EJ  5 

-0.246E-6 

EJ6 

0.558E-6 

EJ7 

-0.326E-6 

J8 

EJ8 

-0.209E-6 

J9 

EJ9 

-0. 094E-6 

Earth's  tesseral  harmonics  C99 

C22 

2. 4369E-6 

s22 

S22 

-1 .4005E-6 

S, 

C31 

2.01 92E-6 

S3I 

S31 

0. 2278E-6 

C32 

C32 

0. 7783E-6 

S32 

S32 

-0. 7552E-6 

C33 

C33 

0. 7387E-6 

S33 

S33 

I . 4343E-6 

C4, 

C4I 

-0.5175E-6 

S4 1 

S4I 

-0.4814E-6 

C42 

C42 

0. 3444E-6 

*42 

S42 

0.7021E-6 

C43 

C43 

1 .0390E-6 

S43 

S43 

-0. 1 192E-6 

C44 

C44 

-0. 1846E-6 

S44 

S44 

0. 2508E-6 

Earth's  angular  velocity  (rad/s) 

E OMEGA 

7.2921 15147E-5 

3 -2 

Sun's  gravitational  constant  (km  s ) 

SUNMU 

1327127E5 

3-2 

Moon's  gravitational  constant  (km  s ) 

SELMU 

4902.756 

Earth's  precession  rate  (rad/s) 

PREC 

6.079E-12 

6 PROGRAM  DESCRIPTION 

The  program  starts  by  reading  all  the  input  data  and  carrying  out  any 
necessary  transformations.  The  latter  define  the  orbit  as  the  cartesian 
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components  of  the  geocentric  position  and  velocity  vectors  at  epoch;  if  a 
covariance  matrix  is  to  be  propagated,  this  may  also  have  to  be  transformed  so 
that  its  elements  refer  to  position  and  velocity  components.  Certain  flags  are 
set  and  those  constants  which  are  a function  of  the  input  data  are  calculated. 

A simple  heading  giving  the  perturbing  forces  included  (i.e.  air  drag, 
luni-solar  perturbations)  is  written  to  the  line  printer. 

The  integration  period  and  interpolation  interval  are  converted  to  days 
and  the  time  of  the  first  interpolation  found  relative  to  epoch.  The  integration, 
which  may  be  forwards  or  backwards  in  time,  is  begun,  and  after  each  integration 
step  the  required  interpolations  are  performed.  The  time  of  the  first  interpola- 
tion occurring  after  the  latest  step  is  then  set.  After  each  interpolation, 
according  to  the  output  options  chosen,  the  interpolated  position  and  velocity 
components  may  be  written  directly  to  the  line  printer  (and  if  required  to  a 
permanent  disc  file)  or  transformed  into  osculating  elements  and/or  range,  range 
rate,  azimuth  and  elevation  for  a given  ground  station.  These  latter  forms  are 
stored  consecutively  in  an  unformatted  scratch  disc  file.  The  associated 
covariance  matrix  may  also  be  written  directly  to  the  line  printer. 

At  each  integration  step,  the  rate  of  change  of  the  radius  vector  is 
calculated.  This  quantity  changes  sign  at  an  apse  enabling  the  integration 
steps  on  either  side  to  be  identified.  Interpolations  are  performed  to  find  the 
time  of  the  apse  and  the  corresponding  osculating  elements  and  radial  distance. 
This  information  is  then  written  to  a second  unformatted  scratch  disc  file. 

When  the  integration  is  complete,  each  ephemeris  is  read  back,  one  line  at 
a time,  from  the  scratch  files  and  output  on  the  line  printer.  Finally,  the 
date,  time  (MJD  and  fraction  of  a day)  and  components  of  position  and  velocity 
for  the  last  integration  step  are  printed.  This  information  may  be  used  to 
restart  the  integration  at  a future  time.  A sample  line  printer  output  is  given 
in  Fig. 3. 

7 DATA  DECK 

7.  1 Overall  description 

This  section  describes  the  input  data  for  POINT,  which,  for  convenience, 
is  treated  here  as  a set  of  punched  cards.  There  are  24  possible  records  which 
may  anpear  in  a data  deck,  and  the  number  present  for  any  one  run  will  be  not 
less  than  4.  The  order  of  the  cards  is  as  follows: 
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(i) 

time  card 

(ii) 

control  cards 

(iii) 

station  cards 

(iv) 

orbit  information  cards 

(v) 

covariance  matrix. 

Items  (iii)  and  (v)  are  options  and  may  be  omitted  in  appropriate  circumstances. 
All  records  are  read  in  free  format  except  for  those  specifying  RAE  or  SDC 
elements  and  those  specifying  a covariance  matrix.  Listings  of  specimen  data 
decks  are  given  in  Fig. A.  (The  first  deck  shown  was  used  to  generate  the  output 
given  in  Fig. 3.)  In  the  descriptions  of  the  data  cards  which  follow,  all  para- 
meters and  constants  are  referred  to  by  their  Fortran-variable  names. 

7. 2 The  time  card 

The  time  card  contains  five  parameters  specifying  the  epoch  of  the  orbital 
elements,  the  time  span  of  the  integration,  the  interpolation  interval  and  the 
time  at  which  the  first  ephemerides  are  required. 

MJD  This  specifies  the  MJD  number  of  the  epoch  at  which  the  orbital 

parameters  of  the  satellite  are  defined. 

EP  This  is  the  time  of  (the  input)  epoch,  expressed  as  a fraction  of 

a day  relative  to  MJD. 

Both  MJD  and  EP  are  set  to  zero  if  the  orbit  is  specified  by  SDC  elements  since 
the  epoch  is  specified  in  the  standard  format  (read  by  subroutine  SDC2EL) . 

TIOR  This  specifies  the  time-interval  of  the  epheraeris,  in  hours. 

DI  This  is  the  required  interpolation  interval,  in  minutes  (negative  if 

the  integration  is  backwards  in  time). 

TINT  This  is  the  time  at  which  the  first  ephemerides  are  required  in 

hours  from  epoch.  If  TINT  = 0.0  , it  is  assumed  to  be  DT  minutes 
from  epoch. 

7. 3 First  control  card 

This  control  card  contains  eight  parameters  (all  integer),  which  control 
the  working  of  the  program. 

I This  indicates  the  form  in  which  the  orbital  parameters  are  provided 

(see  section  7.6).  If  I = 1 , a set  of  injection  conditions  in 
kilometres,  kilometres  per  second  and  degrees  is  read  in  the  order 
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specified  in  section  2.  If  I = 2 , a set  of  position  and  velocity 
components  in  kilometres,  kilometres  per  second  is  read.  If  1=3 
a set  of  osculating  elements  in  kilometres  and  degrees  is  provided. 

If  1=4,  SDC  two-line  elements  are  read.  If  1=5,  RAE  mean 
elements  are  provided. 

NM  This  indicates  the  matrix  operations  to  be  performed  (see  section  7.7). 

If  NM  = 0 , there  are  no  matrix  operations.  If  NM  = 1 , a covariance 
matrix  in  terms  of  injection  conditions  (with  units  feet,  feet  per 
second,  degrees)  is  read,  and  a matrix  in  terms  of  position  and 
velocity  components  (in  units  kilometres,  kilometres  per  second) 
derived  from  it.  If  NM  = 2 , a covariance  matrix  in  terms  of  posi- 
tion and  velocity  components  in  kilometres,  kilometres  per  second  is 
read . 

NOTPT  If  NOTPT  > 0 , an  ephemeris  of  time  and  corresponding  position  and 

velocity  components  (in  days,  kilometres,  kilometres  per  second)  is 
written  to  an  unformatted  disc  file  at  constant  time  intervals  (DT) . 

NPV  If  NPV  = 1 , an  ephemeris  of  time,  position  components,  radial  dis- 

tance, velocity  components  and  speed  is  written  to  the  line  printer. 

NEL  If  NEL  = 1 , an  ephemeris  of  time,  osculating  elements  and  radial 

distance  is  written  to  the  line  printer. 

NAP  If  NAP  = 1 , a table  of  apses  giving  time,  osculating  elements  and 

radial  distance  is  written  to  the  line  printer. 

NAE  If  I i NAE  £ 4 , tables  of  look  angles  for  NAE  stations  giving 

times,  azimuth,  elevation,  range  and  range  rate  are  written  to  the 
line  printer. 

LSP  If  LSP  = I , luni-solar  perturbations  are  included  in  the 

integration. 

7.4  Second  control  card 

This  card  contains  one  control  parameter,  one  parameter  which  is  used  both 

as  a control  parameter  and  a data  item  and  two  data  parameters. 

CJ[  Although  the  integration  step  length  (H)  is  set  in  the  BLOCK  DATA 

segment,  it  may  be  modified  to  H/Cl  by  setting  Cl  #0  (Cl  must  be 
negative  if  backward  integration  is  required). 
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1 1 

2 

ARMA  This  is  the  area-to-mass  ratio  of  the  satellite  (m  /kg).  If  set  to 

0.0,  the  air  drag  terms  are  excluded  from  the  integration. 

SMEAN  This  is  the  value  of  the  solar  10.7cm  radiation  flux  averaged  over 

-22  -2  - 1 

3  days  and  measured  in  units  of  10  W m Hz  . If  ARMA  = 0.0  , 
then  SMEAN  = 0.0  . 

SOLBAR  This  is  the  value  of  the  solar  10.7cm  radiation  flux  averaged  over 

. . -22  -2  -1 

3 months  and  measured  in  units  of  10  W m Hz  If  ARMA  = 0.0  , 

then  SOLBAR  = 0.0  . 

7. 5 Station  cards 

Tables  of  look  angles  may  be  produced  for  NAE  stations,  where  NAE  < 4 . 
One  card  is  provided  for  each  station  and  has  four  parameters. 

SLT  This  is  the  station's  geodetic  latitude  (degrees). 

SLO  This  is  the  station's  geocentric  longitude  (degrees). 

R This  is  the  station's  height  (m)  above  the  geoid. 

STANAM  This  is  the  station  name  of  no  more  than  20  characters. 

7.6  Orbit  information  cards 

Orbit  information  may  be  supplied  in  one  of  five  forms,  each  of  which 
requires  one  or  more  cards.  The  form  is  indicated  by  the  parameter  I on  the 
first  control  card. 

I Number  of  cards  Contents 

1 1 Injection  parameters  (km,  km/s,  deg) 

2 1 Geocentric  position  and  velocity  com- 

ponents (km,  km/s) 

3 1 Osculating  elements  (km,  deg) 

4 2 SDC  two-line  elements 

5 5 RAE  mean  elements 

7.7  Covariance  matrix 

The  covariance  matrix  is  supplied  or  twelve  cards.  The  form  is  indicated 
by  the  parameter  NM  on  the  first  control  card.  Each  card  contains  three 
elements  punched  in  3E15.8  format.  Cards  one  and  two  hold  the  first  row,  cards 
.three  and  four  the  second  row  and  so  on. 
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8 DESCRIPTION  OF  OUTPUT 

A specimen  of  the  line  printer  output  is  shown  in  Fig. 3.  Time  is  given  in 
days  and  fractions  of  a day,  distances  in  kilometres  and  angles  in  degrees. 

A total  of  six  output  options  are  available: 

(1)  An  ephemeris  of  time,  position  components,  radial  distance,  velocity 
components  and  speed  written  to  the  line  printer. 

(2)  A propagated  covariance  matrix,  printed  after  each  interpolation. 

(3)  An  ephemeris  of  time,  position  and  velocity  components  written  to  an 
unformatted  disc  file. 

(4)  An  ephemeris  of  time,  osculating  elements  and  radial  distance  written  to 
the  line  printer. 

(5)  Tables  of  look  angles  for  a maximum  of  four  stations.  Each  table  comprises 
time,  azimuth,  elevation,  range  and  range  rate;  this  information  is  only 
printed  if  the  elevation  is  greater  than  or  equal  to  0.0°. 

(6)  A table  of  apses  written  to  the  line  printer.  Each  entry  gives  the  time, 
osculating  elements  and  radial  distance  of  the  apse.  If  no  apses  are 
found,  an  appropriate  comment  is  printed. 

Finally  at  the  end  of  each  output,  the  time  of  the  last  integration  step  is 
printed  together  with  the  position  and  velocity  components  at  that  step. 

All  or  any  combination  of  the  options  may  be  chosen.  However,  option  2 is 
meant  to  be  used  with  option  1 since  the  program  is  designed  to  print  the  propa- 
gated matrix  after  each  line  of  option  I . If  option  2 is  chosen  without  option  1 , 
the  program  will  function  correctly,  but  time  will  not  be  printed  with  each 
matrix. 

9 PROGRAMS  DERIVED  FROM  POINT 

POINT  has  been  designed  as  a versatile  program  with  a number  of  input  and 
output  options.  However,  for  some  tasks  which  require  much  computation,  these 
options  may  not  be  required  and  it  may  well  be  more  efficient  in  terms  of  execu- 
tion time  and  core  store  to  produce  a special  version  of  the  program.  The 
on-line  editor  makes  this  a straightforward  task. 

Such  a version  was  produced  for  the  studies  of  the  proposed  Astronomical 
Roentgen  Observatory  (ASRO) 1 . Other  versions  which  currently  exist  are  P0INT2, 

PT ATMAN  and  PTDEC. 
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P0INT2  is  used  to  generate  random  transfer  orbits  for  the  synchronous 

. 1 2 
mission  analysis  program  SYNMAP  and  has  been  described  elsewhere 

PTATMAN  comprises  the  standard  version  of  POINT  plus  the  facility  to 
include  a gross  attitude  manoeuvre,  such  as  would  be  carried  out  in  a transfer 
orbit, to  realign  a spacecraft  for  the  firing  of  an  apogee  boost  motor.  This 
manoeuvre  is  assumed  to  be  impulsive  (and  applied  at  apogee)  and  the  program  is 
supplied  with  transverse,  normal  and  radial  delta-velocity  components.  An 
extra  card,  punched  in  free  format,  is  inserted  after  the  second  control  card  in 
the  standard  POINT  data  deck.  This  card  contains  the  apogee  number  (punched  as 
an  integer)  and  the  three  delta-velocity  components.  One  other  feature  of  this 
version  is  that  the  number  of  stations  for  which  look  angles  may  be  produced  has 
been  increased  from  four  to  five,  since  this  is  the  number  usually  associated 
with  the  launch  phase  of  a synchronous  orbit  mission. 

PTDEC  is  used  to  find  the  time  at  which  a satellite  decays,  i.e.  the  time 
at  which  it  falls  below  a given  altitude.  The  program  is  provided  with  the 
area-mass  ratio  of  the  satellite  and  either  a set  of  geocentric  position  and 
velocity  components  or  a set  of  SDC  two-line  elements.  Air  drag,  luni-solar 
perturbations  and  zonal  harmonics  up  to  are  automatically  included  but 

tesseral  harmonics  are  omitted.  The  orbit  is  integrated  forward,  and  at  each 
perigee,  the  osculating  elements,  perigee  radius  and  preceding  apogee  radius  are 
printed.  The  integration  ceases  when  either  the  satellite's  altitude  falls  below 
a given  value  or  a specified  time  limit  has  elapsed. 
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Appendix  A 

POINT  PROGRAM  UNITS 

reduces  an  angle  to  the  range  0 to  2v 
reduces  an  angle  to  the  range  -it  to  it 
sets  certain  constants 

converts  coordinates  to  osculating  elements 
integrates  differential  equations 

determines  eccentric  anomaly  from  Kepler's  equation 
converts  osculating  elements  to  coordinates 
finds  a named  area  on  a disc  file 

converts  a number  into  its  integer  and  fractional  parts 

computes  the  position  and  velocity  components  and  covariance 
matrix  by  interpolation 

computes  the  position  and  velocity  components  between  integration 
steps  by  interpolation 

converts  injection  parameters  to  coordinates 

performs  matrix  multiplication 

computes  atmospheric  density 

controls  orbit  integration  and  interpolation 

transforms  position  and  velocity  into  required  form  and  stores 

auxiliary  for  DEQRSPT;  computes  perturbing  accelerations 

computes  position  and  velocity  components  from  RAE  mean  elements 

produces  line  printer  output  from  information  stored  on  discs 

computes  range,  range  rate,  azimuth  and  elevation  for  given  ground 
station 

reads  RAE  5-card  orbital  elements 
reads  input  data 
releases  peripheral  unit 

performs  rotations  about  an  axis  through  a given  angle  at  a given 
angular  velocity 

forms  scalar  product 

reads  and  converts  SDC  2-line  elements  to  RAE  mean  elements 

sets  initial  values  of  partial  derivatives 

computes  short-periodic  perturbations 

reads  sun/moon  coordinates  from  disc  and  interpolates 

performs  interpolation 

transposes  a matrix 

determines  polar  coordinates  from  cartesian  (two-dimensional) 
permits  disc/core  transfers 
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(a)  The  subprograms  INFIND,  INTFRC  and  SCPROD  are  written  in  PLAN  and  are  used 
as  semi-compiled  segments. 

(b)  The  subprograms  RLEASE  and  UTD4 , though  not  part  of  standard  FORTRAN,  are 
provided  automatically  by  the  1900  series  FORTRAN  compilers,  and  are  not 
described  here. 
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Appendix  B 

SUBPROGRAM  SPECIFICATIONS 


Summary 


Language 

Author 

Function  statement 
Input  argument 
X 


FUNCTION  ANGLE 

- The  function  reduces  an  angle  x (in  radians)  to  the 
range  0 * x < 2ir  . 

- 1900  Fortran. 

- Diana  W.  Scott  (April  1969). 

- FUNCTION  ANGLE (X). 


Angle  x in  radians. 


Output  function 
ANGLE 

Use  of  COMMON 
Source  deck 


Local  storage  used 
Subordinate  subprograms  - None 


- Value  of  x + 2nu  , such  that  0 < ANGLE  < 2ir 

- None. 

- 6 cards,  including  1 comment  card  (ICL  code). 

- 1 real  variable. 


Explanation  - The  standard  function  AMOD  is  used  to  give  the 

fractional  part  of  x/2ir  . If  this  is  negative,  2ir  is  added  to  the  result. 
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FUNCTION  ANGL1 


Summary 

- The  function  reduces  an  angle  x to  the  range 

-IT  < X £ Tt  . 

Language 

- 1900  Fortran. 

Author 

- Diana  W.  Scott  (April 

1969). 

Function  statement 

- FUNCTION  ANGLl(X) 

Input  argument 

- 

X 

Angle  x in  radians. 

Output  function 

- 

ANGL1 

The  value  of  x ± 2mr 

such  that  -ir  < ANGL1  < it 

Use  of  COMMON 

- None. 

Source  deck 

- 7 cards,  including  one 

common  card  (ICL  code). 

Local  storage  used  - 1 real  variable. 

Subordinate  subprograms  - None. 

Explanation  - The  standard  function  AMOD  is  used  to  give  the 

fractional  part  of  x/2ir  . If  this  is  greater  than  ir  , 2ir  is  subtracted.  If 
it  is  less  than,  or  equal  to,  -it  , 2tt  is  added. 
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BLOCK  DATA  (POINT) 

Summary  - The  Fortran  BLOCK  DATA  segment  is  used  to  set  initial 

values  for  certain  quantities  stored  in  conmon  blocks 
/PETURB/,  /CINTEG/,  /CONST/  and  /CON/. 

Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (January  1975). 

Data  in  /PETURB/ 


Variable  name 

Value 

Explanation 

HALFCD 

1 . 1 

Product  $C  . 

RATE 

1 . 1 

Ratio  of  angular  velocity  of  the  atmosphere 

to  that  of  the  earth. 

Data  in  /CINTEG/ 

- 

Variable  name 

Value 

Explanation 

H 

2ir/ 96 

Integration  step  length. 

PD  (18) 

All  elements  0.0 

Partial  derivatives  of  position. 

PDVEL  (18) 

All  elements  0.0 

Partial  derivatives  of  velocity. 

Data  in  /CONST/ 

- 

Variable  name 

Value 

Explanation 

EMU 

398601. 3kra3s"2 

Earth's  gravitational  constant. 

EJ2 

1 082. 637E-6 

Earth's  second  zonal  harmonic  . 

EJ3 

-2.5310E-6 

Earth's  third  zonal  harmonic  • 

EJ4 

-1 .6190E-6 

Earth's  fourth  zonal  harmonic  . 

C22 

2.4369E-6 

Tesseral  harmonic  coefficient  C22 

• 

S22 

-1.4005E-6 

Tesseral  harmonic  coefficient  S22 

• 

C33 

0. 7387E-6 

Tesseral  harmonic  coefficient  C^ 

• 

S33 

1 .4343E-6 

Tesseral  harmonic  coefficient  S^ 

• 

C44 

-0. 1846E-6 

Tesseral  harmonic  coefficient  C^ 

• 

S44 

0. 2508E-6 

Tesseral  harmonic  coefficient  S^ 

• 

C3I 

2.0I92E-6 

Tesseral  harmonic  coefficient 

• 

S31 

0. 2278E-6 

Tesseral  harmonic  coefficient  S^ 

• 

C42 

0. 3444E-6 

Tesseral  harmonic  coefficient 

• 

S42 

0.7021E-6 

Tesseral  harmonic  coefficient 

• 

ERAD 

6378.163km 

Earth's  mean  equatorial  radius. 
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Variable  name 

Value 

Explanation 

EOMEGA 

7.2921 15147E-5rad/s 

Earth's  angular  velocity. 

PREC 

6.079E-12rad/s 

Earth's  precession  rate. 

SUNMU 

1 . 327 1 27E 1 lkm3s_2 

Sun's  gravitational  constant. 

SELMU 

4902.756kro3s“2 

Moon's  gravitational  constant. 

Data  in  /CON/ 

- 

Variable  name 

Value 

Explanation 

EJ5 

-0.246E-6 

Earth’s  fifth  zonal  harmonic  J,.  . 

EJ6 

0.558E-6 

Earth's  sixth  zonal  harmonic  Jg  . 

EJ7 

-0.326E-6 

Earth's  seventh  zonal  harmonic  J^. 

EJ8 

-0.209E-6 

Earth's  eighth  zonal  harmonic  Jg . 

EJ9 

-0.094E-6 

Earth's  ninth  zonal  harmonic  Jg  . 

C43 

1 .0390E-6 

Tesseral  harmonic  coefficient 

S43 

-0.11 92E-6 

Tesseral  harmonic  coefficient  S._* 

4 J 

C41 

-0.5175E-6 

Tesseral  harmonic  coefficient  C,,. 

4 1 

S41 

-0.48I4E-6 

Tesseral  harmonic  coefficient  S,,* 

4 1 

C32 

0.  7783E-6 

Tesseral  harmonic  coefficient  C.^* 

S32 

-0.7552E-6 

Tesseral  harmonic  coefficient  S^’ 

Source  deck  - 18 

cards  (ICL  code) . 

129 


Appendix  B 


21 


SUBROUTINE  COTOE L 


Summary 

Language 

Author 

Subroutine  statement 

Input  arguments 
X,  Y,  Z 

XDOT,  YDOT, 
ZDCT 

EMU 

Output  arguments 
A 
E 

ORB INC 

RANODE 

ARGPER 

EM 

EN 


- The  subroutine  derives  the  standard  elliptic  (osculat- 
ing) orbital  elements  of  an  earth  satellite,  given  its 
position  and  velocity  components. 

- ASA  Fortran  (Standard  Fortran  4). 

- R.H.  Gooding  (May  19.76). 

- SUBROUTINE  COTOEL  (X,  Y,  Z,  XDOT,  YDOT,  ZDOT,  EMU,  A, 
E,  ORBINC,  RANODE,  ARGPER,  EM,  EN) . 


Geocentric  cartesian  components,  (x,y,z),  of  the  satellite, 
x towards  the  vernal  equinox  and  z towards  the  north 
pole. 

Velocity  of  the  satellite,  (x,y,z),  measured  in  the 
same  coordinate  system. 

Earth's  gravitational  constant,  p . 


Semi -major  axis,  a , in  the  same  units  as  x,  y,  z 
Eccentricity,  e . 

Orbital  inclination,  i . 

Right  ascension  of  the  ascending  node,  ft  . 

Argument  of  perigee,  w . 

Mean  anomaly,  M . 

Mean  motion,  n . 


Use  of  COMMON 


Source  deck 


- None. 

- 30  cards,  including  four  conment  cards  (ICL  code) 

- 13  real  variables. 


Local  storage  used 
Subordinate  subprograms  - The  subroutine  TRINV. 

Explanation  - Although  the  subroutine  essentially  transforms  from  the 


six  components  of  the  position  and  velocity  of  a satellite  to  its  six  orbital 
elements,  a seventh  input  argument  permits  an  arbitrary  value  of  the  constant  p 
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to  be  used;  the  seventh  output  argument,  n , is  derived  from  a and  y by  the 
2 3 

relation  n a = y (Kepler's  third  law).  This  means  that  the  subroutine  may  be 
used  very  generally;  e.g.  for  a planet,  taking  the  sun's  y and  interpreting  ft 
as  celestial  longitude.  Units  of  time  and  distance  are  arbitrary  but  must,  of 
course,  be  consistent;  all  angles  are  in  radians. 


The  subroutine  has  been  written  to  give  maximum  accuracy.  All  angles  are 

derived  from  knowledge  of  both  sine  and  cosine,  and  in  such  an  order  that  there 

is  no  difficulty  near  the  singularities  at  e = 0,  i = 0 and  i = ir  , The 
• 1 3 

general  method  is  that  of  Brouwer  and  Clemence  (section  27  of  chapter  1). 


The  first  quantities  to  be  derived  are  a,  e and  the  eccentric  anomaly, 

E . (NB  The  Fortran  variable  E refers  to  the  eccentricity  and  not  the  eccentric 
anomaly.)  These  come  from  the  relations 


and 


e cos  E 


1 


e sin  E 


(xx 


(pa)* 


Z7) 


* 


where  r = (x  + 


2 2,4 

y + z )* 


and 


v2  = 


.2 

+ y + 


Then  M follows  at  once,  since 

M = E - e sin  E . 


Note  that  if  e is  close  to  zero  E is  ill-determined,  reflecting  the 
indeterminancy  of  perigee  - and  in  fact  if  e = 0,  E is  set  to  0 . This  does 
not  matter  at  all,  since  whatever  position  of  perigee  is  specified  by  the  value 
taken  for  E , the  value  of  M is  fully  consistent  with  it. 

The  values  of  i,  ft  and  u>  are  derived  from  the  basic  matrix  relation 
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where  P = (x/r)  cos  E - x(a/y)2  sin  E , 
x 

Q ■ (1  - e2)  ^ [(x/r)  sin  E + x(a/u)^(cos  E - 

Rx  = (ua(l  - e2))  ^ (yz  - zy)  , 

and  similarly  for  P , P , Q„ , Q , R and  R 
J y z y z’  y z 

Thus  i and  ft  are  derived  from 


sin  i sin  ft  = R , 
x 

- sin  i cos  ft  = R , 

y 

and 

cos  i = R , 
z 

with  an  indeterminacy  in  ft  when  i is  close  to  0 or  it  (ft  is  set  to  0 if 
sin  i = 0). 

Finally  u is  determined  from 

sin  i cos  w = Q 

z 

and 

sin  i sin  u>  = P 


where  we  have  to  be  sure  that  no  difficulty  arises  over  the  e-singularity  or 
either  of  the  i-s ingulari ties . 

There  is  no  trouble  near  the  e-singularity  since  the  quantities  e cos  E 

and  e sin  E are  used  (through  Q and  P ) , in  the  derivation  of  w , in  the 

z z 

same  ratio  as  in  the  derivation  of  E , so  that  E + w is  always  correct. 

(When  E is  set  conventionally  to  0 because  e = 0 , e cos  E is  set  to  1 (!) 

tc  ensure  the  correct  ratio  between  Q and  P .) 

z z 

In  the  same  way,  w has  to  be  compatible  with  ft  near  the  i-singularity , 
and  this  is  automatic  when  z and  z are  not  both  exactly  zero,  since  ft  and 
a)  both  depend  on  the  ratio  of  these  two  quantities.  When  z and  z cere  both 
exactly  zero,  the  correct  value  of  u>  is  achieved  by  replacing  P^  and  by 

P and  Q , with  an  additional  factor  to  ensure  that  w is  in  the  correct 

y y 

quadrant . 
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SUBROUTINE  DEQRSPT 

Summary  - The  subroutine  integrates  a set  of  up  to  22  simultaneous 

second-order  differential  equations  of  the  form 

y i = fjU.y, yN,  y, yN) , i = , using 

a Gaussian  eighth-order  second-sum  predictor-corrector'^ 
The  integration  is  started  using  a Butcher  sixth-order 
seven  stage  Runge-Kutta  process'^  to  set  up  the  differ- 
ence table  required  before  the  second-sum  procedure  can 
take  over. 


Language 

Authors 

Subroutine  statement 


- 1900  Fortran. 

- A.W.  Odell  and  G.J.  Davison  (April  1973). 

- SUBROUTINE  DEQRSPT  (NTRY,  DAUX) . 


Input  arguments 
NTRY 


1 for  an  initialization  entry  to  the  subroutine,  and 

2 for  all  normal  entries  (see  Explanation). 


DAUX 

Output  arguments 
Use  of  COMMON 


The  auxiliary  subroutine  which  evaluates  y^  . 

- None. 

- Certain  quantities  in  common  block  /CINTEG/  are  used 
as  follows :- 


Input  arguments  in  /CINTEG/  - 

NEQ  Number  of  equations  (normally  4 or  22). 

H Integration  step  size,  h (positive  or  negative). 

Input  and  output  arguments  in  /CINTEG/  - 

T Independent  variable,  t . 

Y(22)  Dependent  variables,  y^  . 

YP(22)  First  derivatives,  y^  . 

Y2P(22)  Second  derivatives,  y.  , as  computed  by  the  auxiliary 

subroutine  DAUX. 


Output  arguments  in  /CINTEG/  - 

IAUX  Set  to  -1  following  the  initialisation  entry  (NTRY=1), 

to  0 if  T has  been  changed  (during  a normal  entry), 
and  to  1 otherwise. 
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Source  deck  - 218  cards  (ICL  code). 

Local  storage  used  - 10  integer  variables,  1 logical  variable,  63  real 

variables  and  814  real-array  elements. 

Subordinate  subprograms  - The  auxiliary  subroutine,  named  in  the  call  to 

DEQRSPT  which  takes  the  role  of  DAUX. 

Explanation  - DEQRSPT  integrates  a set  of  N(<22)  simultaneous 

second-order  differential  equations  of  the  form  y.  = f ^ ( t ,y  ( , . . . y^  ,y  ^ , . . . y„^)  , 

i = 1 , . . . ,N  , using  a Gaussian  eighth-order  second-sum  predictor  corrector.  The 
integration  is  started  using  a Butcher  (6,7)  R-K  process  to  set  up  the  difference 
table  required  before  the  second-sum  procedure  takes  over. 


Prior  to  any  integration,  DEQRSPT  must  first  be  called  with  NTRY  = 1.  This 
is  the  initialization  entry  in  which  the  step  size  h'(=  |h)  is  set  for  the 
Butcher  integration  and  DAUX  is  called  to  evaluate  y\  at  t^  . All  subsequent 
entries  are  made  with  NTRY  = 2,  and  one  integration  step  (two,  whilst  still  in 
the  Butcher  mode)  is  performed  before  control  is  returned  to  the  calling  program. 


If  y.  , and  y.  , are  the  values  of  y.  and  y.  at  t = t the 

lfK  1 * K,  1 1 K 

formulae  used  in  the  next  Butcher  step  to  evaluate  y.  , v.  and 

T.  i , k+  I i , k+ 1 

y.  ...  , the  values  of  y. , y.  and  y.  at  t = t + h'  , are: 

* i " I 1 K 


ri,k+l 


\ \ 

y.  , + ) w k. 

'i.k  s it 


s-  I 


/ 

yi,k+i  - ^i.k  + Yj  Wa*i£ 


s- 1 


and 


fi,k-H 


- f 


i(tk + ^.k+i’yi.k+i’"-^.^  • 1-1 N 


where  k.  " h'f.  t.  ♦ c h',y.  + ) a I , 

is  l k s l ,k  sj  ij  is  1 
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and  the  remaining  quantities  are  given  by  the  following  table 


Given  y^  Q and  y^  q , y^  q is  obtained  from  DAUX,  and  the  first  two 

Butcher  steps  of  size  h'  then  give  y.  , y.  and  y.  at  time 

1,1  1,1  # i,l 

li  " tn  ^ h * The  Proce8s  is  repeated  until  y.  D,  y.  0 and  y.  „ are 

1,01,0  1,0 

obtained  after  a total  of  16  steps  (eight  calls  to  the  subroutine).  The 
Butcher  process  is  now  complete. 

On  the  ninth  call  (with  NTRY  - 2)  to  DEQRSPT  the  following  difference  table 
is  constructed  for  each  y^  = f^  . 
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7? 

1,8  ^ 


i,9 


First  the  known  y.  are  differenced  to  give  the  part  of  the  table  to  the  right 

••  1 . -1-2 
of  the  y , and  then  the  first  and  second  sums  7.  . and  7.  _ are  formed, 

i,4  1,3 

using  the  equations 


7?2,  - h"2y.  . - B.y.  . - B_7?  , - B.7*  , - B,V?  . - B 

1,3  7 1,4  07 1 ,4  2 i,5  4 i,6  6i,7  a 1,8 
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and 

V7\  - h_,y.  . - Any.  , - A,?!  - A.V?  - A,V?  , - A.V*  , 

i,4  7i,4  07i,4  1 i,5  2 i,5  3 i,6  4 i,6 

- A..V5  , - A.V6.  _ - A,V?  „ - A.V?  . . 

5 i,7  6 i,7  7 i,8  8 i,8 

-2  -1  . . 

The  remaining  V and  V quantities  above  the  dotted  line  are  defined  on 
the  basis  that  the  difference  between  any  entry  and  the  entry  above  must  equal 
the  entry  on  the  right, 
and  these  are  given  by 


. -1  -2 

the  entry  on  the  right.  Actually  only  V.  and  V.  are  required  explicitly 

1 , o 1 , o 


vi|8  • v 


lU  * *1,5  * *1.6  * *1.7  * *i,B 


and 


v7 2 - v72.  + 5V71.  + 4y,  , + 3y . , + 2y.  _ + y. 

1,8  1,3  1,4  i,5  7i  ,6  7i  ,7  7i,- 


8 ' 


A table  containing  the  coefficients  used  in  the  above  and  following  equations  is 
appended. 

The  ninth  normal  call  to  the  subroutine  continues  with  an  integration  step, 
using  the  Gaussian  (predictor)  formulae:- 


and 


yi,9  “ h 


+ C'y.  _ 

♦ c.v.  n + 

2 

CV7  „ + 

3 

c„  v:  „ 

07 1 ,8 

1 i,8 

2 i ,8 

3 i,8 

+ c,v?  _ 

+ c,v^  0 + 

C„V?  „ + 

C0V® 

5 i,8 

6 i ,8 

7 i,8 

8 i,8 

yi,9  “ hl 


+ F.v!  „ + F,V?  D ♦ F„V?  „ + F,V? 


‘(7i!s  + Vi, 8 ' 1 ' i , 8 T r2’i,8  T r3’i,8  T r4vi,8 

+ Vi, 8 + V{,8  * Vl,8  + Fi,s)  • 

which  require  only  the  quantities  immediately  above  the  dotted  line  in  the  table. 


y.  q is  then  obtained  using  DAUX  and  the  row  of  differences  under  the 
diagonal  line  found.  This  row  is  then  used  to  obtain  corrected  values  of  y. 
and  y^  ^ using  the  equations 


i,9 
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»i,9  ' * Vi,9  * Dl7!,9  * D27b  * Vi,9  * Vi,' 

* Vi, 9 * Vi, 9 * Vi,.  * Vb) 


• / -|  1 2 3 4 

y.  0 = h V.  q + E_y . Q + E.V,  Q + E07.  . + E-V7  Q + E.V7  Q 

7 l , 9 \ i , 8 0 l , 9 1 i,9  2 i,9  3 i , 9 4 t , 9 


Vi, 9 * Vi,9  * Vb  * vb)  • 


••  8 

y.  . is  then  redetermined,  and  the  row  of  differences  out  to  V? 

t.9  i, 

the  dotted  line  recalculated.  This  row  is  then  used  to  obtain  the  final 
corrected  values  of  y^  g and  y^  g using  the  above  equations. 

Coefficients  for  A..  B..  C. . D. . E.  and  F. 


5 

6 

191 

191 

60480 

1 20960 

0 

289 

3628800 

275 

33953 

4032 

518400 

19 

9829 

6048 

3628800 

863 

275 

60480 

24192 

19087 

5257 

60480 

17280 

2497 


3628800 


8183 


33953 


3628800 


1070017 


3628800 


2 49  7 


22809600 


3250433 


53222400 


330157 


159667200 


8 


1036800 
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FUNCTION  EAFKEP 


Summary 


Language 

Author 

Function  statement 
Input  arguments 
EM 
ECC 

Output  function 
EAFKEP 


- The  function  solves  Kepler’s  equation;  i.e.  it  provides 
the  eccentric  anomaly  of  a celestial  body,  E , given 
the  orbital  eccentricity,  e , and  mean  anomaly,  M . 
Kepler's  equation  is  M = E - e sin  E . 

- ASA  Fortran  (Standard  Fortran  4). 

- R.H.  Gooding  (May  1976). 

- FUNCTION  EAFKEP  (EM,  ECC). 

Mean  anomaly,  M (radians). 

Eccentricity,  e . 


Eccentric  anomaly,  E . 


Use  of  COMMON 


Source  deck 


- None. 

- 13  cards,  including  3 comment  cards  (ICL  code). 

- 3 real  variables. 


Local  storage  used 
Subordinate  subprograms  - None 
Explanation 

Improved  approximations  are  given  by  Newton's  method;  thus 


- A first  approximation  to  E is  given  by  E^  = M . 


i+ 1 


= E.  + 

l 


(M 


+ e sin  E 


i>/(1  - 


cos  E^) 


The  process  ends  after  three  iterations  if  e < 0.003  , and  otherwise  after  five. 
This  ensures  sufficient  accuracy  at  all  times,  while  making  the  subroutine 
independent  of  the  word-length  of  the  computer.  (If  the  magnitude  of  |e^+)  “ E^| 
were  used  as  a criterion  for  convergence,  the  numerical  value  it  was  compared 
with  would  have  to  vary  from  computer  to  computer.) 
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Summary 

Language 

Author 

Subroutine  statement 

Input  arguments 
A 
E 

El 

BO 

SO 

EM 

EMU 

Output  arguments 
X.Y.Z 

XDOT , Y DOT , ZDOT 

Use  of  COMMON 

Source  deck 

Local  storage  used 

Subordinate  subprograms 

Explanation 


SUBROUTINE  ELTOCO 

- The  subroutine  converts  osculating  Kepler  elements 
into  position  and  velocity  components. 

- ASA  Fortran  (Standard  Fortran  4) . 

- A.W.  Odell  (August  1970). 

- SUBROUTINE  ELTOCO  (A,  E,  El,  BO,  SO,  EM,  EMU,  X,  Y,  Z, 
XDOT,  YDOT , ZDOT). 

Semi-major  axis,  a . 

Eccentricity,  e . 

Orbital  inclination,  i . 

Right  ascension  of  ascending  node,  12  . 

Argument  of  perigee,  u>  . 

Mean  anomaly  M if  EMU  'i  0 , and  true  anomaly  v 
if  EMU  < 0 . 

Earth's  gravitational  constant,  u , the  absolute 
value  is  used. 


The  geocentric  cartesian  position  coordinates  (x,y,z) 
of  the  satellite,  x towards  the  vernal  equinox  and 
z towards  the  north  pole. 

Velocity  coordinates  (x,y,z)  of  the  satellite  measured 
in  the  same  coordinate  system. 

- None. 

- 26  cards  (ICL  code). 

- 7 real  variables. 

- The  function  EAFKEP  and  the  subroutine  ROTATE. 

-If  M is  supplied,  it  is  converted  to  the  eccentric 


anomaly  E by  the  function  EAFKEP,  and  sin  v and  cos  v are  computed  using: 
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2 j / 

sin  v = (1  - e )J  sin  E / ( 1 - e cos  E) 

and 

cos  v = (cos  e - e)/(l  - e cos  E)  . 

Given  sin  v and  cos  v , r,  x,  y,  z,  x,  y,  z are  computed  using: 

r = a(l  - e2)/(l  + e cos  v) 
x = r cos  v , 

y = r sin  v , 

z = 0 , 

x = - [p/(a(I  - e2))]^  sin  v , 
y = [u/(a(l  - e ))]*  (e  + cos  v) 

and 

z = 0 . 


Three  rotations  are  then  performed  using  subroutine  ROTATE,  one  about  the  normal 
to  the  orbital  plane  (z-axis)  by  an  angle  io  , one  about  the  nodal  line  (x-axis) 
by  an  angle  i and  one  about  the  polar  axis  (z-axis)  by  an  angle  ft  , giving 
the  required  result. 

If  M is  supplied  then  e should  be  <1.  If  v is  supplied  this  limita- 
tion does  not  apply. 
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Summary 

Language 

Author 

Function  statement 
Input  argument 
NAME 


FUNCTION  INFIND 

- The  function  finds  the  location  of  a named  area  on  a 
specified  disc  file. 

- PLAN  for  use  with  1900  Fortran. 

- A. W.  Odell  (July  1973). 

- FUNCTION  INFIND  (NAME). 


Name  of  an  area  on  the  disc  file  up  to  12  characters  in 
length;  must  be  array  element  or  text  (Hollerith) 
constant. 


Output  function 
INFIND 


Number  specifying  an  area  on  the  disc  file,  0 if  the 
name  is  not  found  in  the  index. 


Use  of  COMMON  - The  first  128  integer  locations  of  blank  common  are 

used  as  temporary  working  space. 

Source  deck  - 37  cards,  including  3 comment  cards  (IBM  'bed'  code). 

Local  storage  used  - 37  words  for  program,  7 words  for  data. 

Subordinate  subprograms  - The  subroutine  UTD4. 

Explanation  - The  subroutine  assumes  that  an  index  has  been  set  up  on 

the  disc  file  using  subroutine  INITD  and  that  information  has  been  put  in  the 
index  using  subroutine  ADDINF.  Associated  with  each  name  in  the  index  is  a 
number  specifying  an  area  on  the  disc  file.  If  this  name  is  not  found,  INFIND 
is  set  to  0;  otherwise  it  is  set  to  the  associated  number. 

Remark:  A Fortran  version  of  this  subroutine  exists. 
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FUNCTION  1NTFRC 


Summary 

- The  function  computes  the  integral 
of  a number. 

and  fractional 

parts 

Language 

- PT.AN,  for  use  with  1900  Fortran. 

Author 

- A.W.  Odell  (February  1971). 

Function  statement 

- FUNCTION  INTFRC  (X) . 

Input  and  output  argumen 

ts  - 

X 

Real  number,  which  is  truncated  to 

its  fractional 

part . 

Output  function 

- 

IN’TFRC 

Integral  part  of  X. 

Use  of  COMMON 

- None. 

Source  deck 

- 23  cards  including  3 comment  cards 

(ICL  code). 

Local  storage  used 

- 18  words  for  program. 

Subordinate  subprograms 

- None. 

Exp  1 anation 

fractional  part.  For  ex 

- X is  split  into  its  mathematical  integral  part  and 
ample : 

its 

X = - 3.4 

I = INTFRC  (X) 

would  result  in  I being  set  to  -4  and  X to  0.6. 


A Fortran  version  of  this  subroutine  is  also  available. 
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Summary 


Language 

Author 

Subroutine  statement 
Input  arguments 
TP 

TIM 

PPVA(9) 

PD1 (18) 

PDV()8) 

PDA (1 8) 

ABSTIM 
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SUBROUTINE  INTPTA 

- The  subroutine  calculates  the  cartesian  components  of 
the  geocentric  position  and  velocity  of  an  earth  satel- 
lite at  constant  time  intervals  by  interpolation 
between  two  given  points  in  the  orbit.  These  points 
are  defined  by  position,  velocity  and  acceleration 
components  at  specified  times.  If  required,  the 
subroutine  will  also  interpolate  for  the  elements  of 
the  covariance  matrix. 

- 1900  Fortran. 

- M.D.  Palmer  (July  1974) 

- SUBROUTINE  INTPTA  (TP,  TIM,  PPVA,  PD1 , PDV,  T2,  ABSTIM). 


Time  (tj)  prior  to  the  first  interpolation,  measured  in 
days  from  epoch. 

Time  (t^)  after  the  first  interpolation,  measured  in 
days  from  epoch. 

(When  integrating  backwards,  TP  and  TIM  will  both 
be  negative.) 


Satellite  geocentric  position,  velocity  and  accelera- 
tion components  at  time  t^ 
(x1,y1,z],x1,y1,z1,X),y),21)  = (rj.Vj.v,)  . 

Partial  derivatives  of  position  at  time  tf  , 

Partial  derivatives  of  velocity  at  time  t^  , 

3 V3  (vy 

Partial  derivatives  of  acceleration  at  time  t| 

8Vi/8(ro,v0)  . 


The  absolute  value  of  t ^ . 
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Output  arguments 
T2 

Use  of  COMMON 


The  absolute  value  of  the  time  of  the  first  interpola- 
tion after  t^  . 

- Certain  quantities  in  common  blocks  /CINTEG/  and 
/INTP/  are  used  as  follows: 


Input  arguments  in  /CINTEG/  - 


MJDOCH 
N ( I ) 

PP  (3) 

PD (18) 
PDVEL( 18) 


Modified  Julian  day  number  of  epoch. 

N(l)  =4  if  only  the  position  and  velocity  components 
are  required  and  N(l)  = 22  if  the  covariance  matrix 
is  also  required. 

Satellite's  geocentric  position  components  at  time  t9 
r2  — (X2*^2,^2^" 

Partial  derivatives  of  position  at  time  - 

3f2/3(l£’V  * 

Derivatives  with  respect  to  the  independent  variable 
s of  the  partial  derivatives  of  position  at  t9  - 


ds 


3 r- 


3(ro*vo) 


TVEL 
PV  (3) 

PA(3) 

PDACCT ( 1 8) 


Value  of  dt/ds  at  time  t^  . 
Satellite's  velocity  components  at  time 


t2  9 v2 


(x2 >y2> z2^ * 


Satellite's  acceleration  components  at  time 
t 2 t v 2 - (x2>y2»z2) • 

Partial  derivatives  of  acceleration  at  time  t. 


3v2/3(r0,vo)  • 


Input  arguments  in  /INTP/  - 


EP 


FDT 


Time  of  epoch,  expressed  as  a fraction  of  a day  relative 
to  MJDOCH. 

Interval  between  interpolations  (days). 

Covariance  matrix  in  terms  of  the  geocentric  position 
and  velocity  components  at  epoch,  cov  (r^.v^)  . 


COVX 
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NOTPT  If  NOTPT  > 0 , the  interpolated  output  is  stored  in  a 

disc  file. 

Input  and  output  arguments  in  /INTP/  - 

TINT  On  input,  the  time  (t)  of  the  first  interpolation 

measured  in  days  from  epoch.  On  output,  the  time  of 
the  first  interpolation  after  t^  , in  days  from  epoch. 
(If  integrating  backwards,  TINT  is  negative.) 

Subordinate  subprograms  - The  subroutines  COTOEL,  MATMUL,  OUTPUT,  RAZEl,,  TRAMAT 

and  TRINV  and  the  functions  ANGLE  and  INTFRC. 


Local  storage  used  - 108  real  array  elements,  24  real  variables  and  3 integer 

variables . 


Source  deck 


- 68  cards  (ICL  code). 


Explanation  - The  subroutine  is  designed  to  carry  out  a series  of 

interpolations  at  constant  time  intervals  between  consecutive  steps  in  an  orbit 
integration.  After  the  interpolated  position  and  velocity  components  have  been 
obtained,  OUTPUT  is  called  to  execute  the  output  options  available  with  POINT. 
If  NOTPT  > 0 , the  MJD, fraction  of  a day,  position  and  velocity  components  are 
stored  in  an  unformatted  disc  file  (channel  5). 

The  formulae  used  for  interpolation  are: 


f(t) 


(1  + 3p  + 6p^)f(tj)  + hp ( 1 + 3p)f'(tj)  + P f"(tj) 


(1  + 3q  + 6q  ;f (t„)  - hq(l  + 3q)f'(t,) 


+ q2f"(t2) 


for  position  components,  and 

- 30h-,p2q2[f(t2)  - f(t,)]  + q2 ( 1 + 5p) ( 1 - 3p)f ' (tj) 

+ p2(l  + 5q) ( 1 - 3q)f*(t2) 

+ J pq2(2  - 5p)f " (t, ) - ^ p2q  (2  - 5q)f"(t2) 


f'(t) 
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for  velocity  components,  where  h“t2-t|,p=(t-  t^/h,  q = 1 - p and 
f (t.),  f'(t^),  f"(t.)  are  the  position,  velocity  and  acceleration  components  at 
time  t^  . 

The  subroutine  will,  if  required,  interpolate  for  the  partial  derivatives 
of  position  and  velocity,  and  use  them  to  obtain  a covariance  matrix,  X , at 
the  required  time  by  means  of  the  transformation 

T 

X = A cov(r0,vQ)  A 


where  A 


3x 

Dx 

3x 

dX 

3x 

3x 

9xo 

9yo 

9zo 

3xq 

9yo 

9zo 

3 1- 

3 L. 

iiL. 

3x0 

9yo 

9zo 

5x0 

9yo 

3zo 

3 z 

3z 

3 z 

3 z 

3z 

3 z 

3xo 

9yo 

3zo 

9xo 

9yo 

3x 

3x 

3x 

3x 

3x 

3x 

9x0 

ayo 

9zo 

3xo 

9yo 

9zo 

3 L 

iL 

iL 

ii_ 

3x0 

9yo 

9zo 

3xo 

9y0 

3zo 

2JL 

3 i 

3 2 

32 

3 i 

3z 

_3xo 

9y0 

9zo 

3i0 

9y0 

1 

o 

• N 

<T> 

and  c0V(tQ,VQ)  is  the  covariance  matrix,  in  terms  of  the  radial  distance 
Tq  , and  satellite  velocity,  Vq  , at  the  start  of  the  integration. 

A is  obtained  by  interpolation  between  the  partial  derivatives  of  position 
and  velocity  at  tj  and  t2  , the  computation  requiring  the  partial  derivatives 
of  acceleration  at  those  times.  The  partial  derivatives  of  position  are  held  in 
the  arrays  PD1  and  PD  respectively,  the  order  of  storage  being 


J 
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PD( 1 )-PD(6)  9x/3xq,  3x/3yQ,  3x/3zQ,  3x/3xQ,  3x/3y0>  3x/3zQ 

PD(7)-PD( 1 2)  3y/3xQ,  3y/3yQ,  3y/3zQ,  3y/3xQ,  3y/3yQ,  3y/3zQ 


PD ( 1 3)-PD( 1 8)  3z/3xq,  3z/3yQ,  3z/3zQ,  3z/3xQ,  3z/3yQ,  3z/3zQ 

and  similarly  for  PD(1). 

The  partial  derivatives  of  velocity  at  tj  are  held  in  the  array  PDV,  the 
order  of  storage  being. 


PDV ( 1 )-PDV (6) 

3x/3Xq  , 

etc . 

PDV ( 7)-PDV ( 1 2) 

3y/3xQ  , 

etc. 

PDV ( 1 3)-PD ( 1 8) 

3 z/3Xq  , 

etc. 

The  partial  derivatives  of  velocity  at  are  not  available  explicitly 

but  are  found  from 


jd_ 

d s 


3r ' 


3(r0’V 


dt 

ds 


where  s is  the  independent  variable  defined  by  dt/ds  = r u * . dt/ds  is 
the  Fortran  variable  TVEL  and 

d_ 

ds 


is  the  array  PDVEL. 

The  partial  derivatives  of  acceleration  at  t;  and  t^  are  held  in  the 
arrays  PDA  and  PDACCT  respectively,  the  order  of  storage  being 

PDA (1  )-PDA (6)  3x/3xq  , etc. 

PDA (7) -PDA ( 1 2)  3y/3xQ  , etc. 

3z/3Xq 


PDA ( 1 3) -PDA (18) 


, etc. 
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and  similarly  for  PDACCT. 

The  matrix  X is  written  to  the  line  printer  and,  if  NOTPT  > 0,  it  is  also 
written  to  the  unformatted  disc  file. 


i 
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Summary 

Language 

Author 

Subroutine  statement 
Input  arguments 
TINT 
TP 

TIMET 

PPVA(9) 

Output  arguments 
PT(3) 

VT(3) 

Use  of  COMMON 

Input  arguments 
PP  (3) 

V(3) 

XACCT(3) 


SUBROUTINE  INTPTB 

- The  subroutine  calculates  the  cartesian  components  of 
the  geocentric  position  and  velocity  of  an  earth  satel- 
lite at  a given  time  by  interpolation  between  two 
points  in  the  orbit.  These  points  are  defined  by 
position,  velocity  and  acceleration  components  at 
specified  times. 

- 1900  Fortran. 

- M.D.  Palmer  (July  1974) 

- SUBROUTINE  INTPTB  (TINT,  TP,  TIMET,  PPVA,  PT,  VT) . 

Time,  t , at  which  the  interpolation  is  required. 

Time,  tj  , of  the  point  prior  to  the  interpolation. 

Time,  t ^ , of  the  point  after  the  interpolation. 

Satellite's  geocentric  position,  velocity  and  accelera- 
tion components  at  time  ^ (Xj,yj,z  etc.). 


Satellite's  geocentric  position  components  at  time  t 
(x,y,z) . 

Satellite's  geocentric  velocity  components  at  time  t 
(x,y,z) . 

- Certain- arguments  in  the  common  block  /CINTEG/  are  used 
as  follows: 


Satellite's  geocentric  position  components  at  time  t-, 
(*2’^  2,Z2^ * 

Satellite's  geocentric  velocity  components  at  time  t0 
(x2,y2,z2) . 

Satellite's  geocentric  acceleration  components  at  time 

C2  ^x2»y2,^2) • 
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Local  storage  used  - 23  real  variables  and  I integer  variable. 

Subordinate  subprograms  - None. 

Source  deck  - 40  cards  (ICL  code) . 

Explanation  - The  subroutine  obtains  a satellite's  geocentric 

position  and  velocity  components  at  a given  time  by  interpolation  between  two 
points  in  the  orbit.  These  will  normally  be  adjacent  steps  in  a numerical 
integration. 


These  equations  used  for  interpolation  are: 


f (t) 


(1  + 3p  + 6p2)f(t  ) + hp(l  + 3p)f ' (t  ) + p2f"(t1) 


(1  + 3q  + 6q  )f(t9)  - hq ( 1 + 3q)f ' (t,) 


+ V cl2f"(t2) 


for  position  components  ai-d 


f'(t)  = 3 Oh  'p2q2[f (t2)  - f(t,)]  + q2 ( 1 + 5p)  ( 1 - 3p)f'(t)) 

+ p2(l  + 5q)  ( 1 - 3q)f'(t2) 

2 

♦ \ ~ 5P)f"(t,)  - $ p2q  (2  - 5q)f"(t2) 

for  velocity  components,  where  h»t2-t)tp-(t-  t()/h,  q - 1 - p and 

f'(t^)  and  f"(t^)  are  the  position,  velocity  and  acceleration  components 
at  time  t.  . 
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SUBROUTINE  1PTOCO 

Summary  - The  subroutine  computes  the  geocentric  cartesian  com- 

ponents of  the  position  and  velocity  of  an  earth  satel- 
lite in  PROP  axes,  given  the  date  and  time  and  a set 
of  the  standard  parameters  used  to  define  orbit 
injection. 

Language  - USA  Standard  Fortran  (USAS  X3.9  - 1966). 

Authors  - G.E.  Cook  and  R.  Clarke  (October  1972). 

Subroutine  statement  - SUBROUTINE  IPTOCO  (X,  Y,  Z,  XDOT , YDOT,  ZDOT,  V,  CA, 

AZ,  R,  XLAT,  XLONG , MJD,  TIME). 


Input  arguments 
V 

CA 

AZ 

R 

XLAT 

XLONG 

MJD 

TIME 


Speed,  v . 

Climb  angle,  6 . 

Azimuth  of  velocity  vector,  T (measured  E of  N) . 
Radial  distance,  r . 

Geocentric  latitude,  <J>  . 

Longitude,  X . 

Modified  Julian  Day  number. 

Time  (fraction  of  a day)  such  that  t is  given  by 
MJD  + TIME  . 


Output  arguments 
X.Y.Z 

XDOT,  YDOT,  ZDOT 
XLONG 

Use  of  COMMON 


Geocentric  coordinates  of  satellite,  x,  v,  z . 
Geocentric  components  of  velocity  vector,  x,  y,  z . 
Longitude  in  inertial  axes,  X'  . 

- None. 

Source  deck  - 22  cards  (ICL  code). 

Local  storage  used  - 8 real  variables. 

Subordinate  subprograms  - None. 

Explanation  - Units  of  length  and  time  are  arbitrary,  except  for  the 


variables  MJD  and  TIME;  angles  are  in  radians. 
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The  components  of  position  and  velocity  are  evaluated  relative  to  the 
following  coordinate  system:  the  origin  0 is  at  the  earth's  centre  and  Oz  points 
towards  the  north  pole;  Ox  lies  in  the  plane  of  the  true  equator  of  date,  but 
instead  of  pointing  towards  the  true  equinox  of  date  it  points  towards  a projec- 
tion of  the  mean  equinox  of  the  epoch  1950.0  (MJD  33281.9234);  Oy  completes  the 
right-handed  system  Oxyz. 

The  injection  conditions  are  specified  relative  to  an  earth-fixed  system 
with  longitude  measured  eastwards  from  the  Greenwich  meridian.  The  longitude  A' 
relative  to  the  inertial  system  defined  above  is  found  by  adding  0 to  A , 

A 

where  0 is  the  'modified  sidereal  angle'  given  by 

0 = 100.075542  + 360.9856 1 2288( t - 33282.0)  , 

i.e.  the  origin  has  been  adjusted  slightly  from  1950.0.  (The  modified  sidereal 
angle  differs  from  sidereal  time  only  because  of  the  choice  of  a non-standard 
reference  direction.) 

The  geocentric  components  of  position  and  velocity  are  obtained  from  the 
following  equations: 

x = r cos  <(>  cos  A' 

y = r cos  sin  A' 

z = r sin 


X » V 

jsin  0 

cos  <j>  cos 

A' 

- cos  0 

(cos 

sin 

$ cos  A + 

sin 

f sin 

\)\ 

y - v 

{sin  0 

cos  4>  sin 

A' 

+ cos  0 

(sin 

cos 

A - cos  V 

sin 

$ sin 

X)\ 

i • v |sin  0 sin  <t>  + cos  0 cos  f cos  $ 
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Summary 

Language 

Author 

Subroutine  statement 
Input  arguments 
EMI 
EM2 
II 
KK 
JJ 
IA 

IB 


SUBROUTINE  MATMUL 

- The  subroutine  performs  matrix  multiplications. 

- USA  Standard  Fortran  (USAS  X3.9  - 1966). 

- G.E.  Cook  (May  1970). 

- SUBROUTINE  MATMUL  (EMI,  EM2 , EM3 , II,  KK,  JJ,  IA,  IB). 

Matrix  of  dimension  (II, KK). 

Matrix  of  dimension  (KK,JJ). 

Actual  number  of  rows  in  EMI  . 

Actual  number  of  columns  in  EMI  and  rows  in  EM2  . 

Actual  number  of  columns  in  EM2  . 

Maximum  number  of  rows  in  EMI  as  given  in  the 
dimension  statement  of  the  calling  segment. 

Maximum  number  of  rows  dimensioned  for  EM2  in  the 
calling  segment. 


Output  arguments 
EM  3 

Use  of  COMMON 

Source  deck 

Local  storage  used 

Subordinate  subprograms 

Explanation 


Matrix  of  dimension  (II, JJ). 

None . 

13  cards,  including  4 comment  cards  (ICL  code). 

3 integer  variables. 

None. 

The  subroutine  performs  the  matrix  multiplication 


M,M2 


= M„ 
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SUBROUTINE  MODAT 


Summary 


- The  subroutine  evaluates  upper-atmosphere  density, 
density  scale  height  and  scale  height  gradient  using 
a simple  analytic  model. 


Language 

Authors 

Subroutine  statement 
Input  arguments 
’IE  IT 
TINF 
TGRADO 


Output  arguments 
DEN 

SCALHT 

SHGRAD 


- 1900  Fortran. 

- G.E.  Cook  and  K.J.  Tomlinson  (April  1969). 

- SUBROUTINE  MODAT  (HE IT,  TINF,  TGRADO,  DEN,  SCALHT,  SHGRAD)  . 

Height  above  the  earth's  surface,  y . 

Exospheric  temperature,  T^  . 

Atmospheric  temperature  gradient  dT/dy  at  the 
reference  altitude,  (120km). 

3 

Atmospheric  density,  p (g/cm  ). 

Density  scale  height  (km). 

Density  scale  height  gradient. 


Use  of  COMMON 
Source  deck 
Local  storage  used 


- None. 

- 46  cards,  including  2 comment  cards  (ICL  code). 

- 24  real  array  elements,  16  real  variables,  1 integer 
variable. 


Subordinate  subprograms  - None. 

Explanation  - The  values  of  density  and  density  scale  height  are 

obtained  from  a simple  analytic  model  of  the  Earth's  upper  atmosphere.  If  gQ 
denotes  the  local  value  of  the  acceleration  due  to  gravity,  the  geopotential 
height  above  the  reference  altitude  y^  is  defined  by 


4 


y 

I |g(y)/g(y0) fdy 


(i) 


so  that 
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c = (y  - y0)(R  + y0>/(R  * y > » (2) 

R being  the  mean  radius  of  the  Earth.  The  temperature  of  the  atmosphere  is 
represented  as  a function  of  geopotential  height  by  the  expression 


T (y)  = T^jl  - a exp(-  t4>| 


(3) 


where 

by 


T 

CO 


is  the  exospheric  temperature  and  a 


a 


T(y0) 

T 


and 


T 


are  constants  defined 


(4) 


and 


x = 


T - T(y 


t(§) 


(5) 


o'  y=yr 


If  the  atmosphere  is  assumed  to  be  in  diffusive  equilibrium  above  the  reference 
altitude,  the  number  density  n^  of  the  ith  constituent  of  molecular  (or  atomic) 
mass  nn  is  given  by 


1 dni 


mi8  1 dT 


n^  dy  kT  T dy 


(1  + a)  , 


(6) 


where  k is  Boltzmann's  constant  and  a is  the  thermal  diffusion  factor. 

With  the  temperature  profile  (3),  equation  (6)  can  be  integrated  to  give 


n. (y«)  ! 


J - a 


ni  = ni(yo;  jr-"a"pFvo 


l+Y.+a 

1 


exp (-  , 


where  Y. 


miR(V 


l xkT 
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The  density  p is  given  by 

P = £ n.m.  . 

i 

For  the  constant  boundary  conditions  at  the  reference  altitude  of  120km 
we  use  the  values  assumed  by  Jacchia*  in  the  construction  of  his  static  diffu- 
sion profiles: 

T ( 1 20)  = 355  K 

n(N  ) = 4.0  x 101 1 cm-3 

n(02)  = 7.5  x io'°  cm"3 

n(0)  = 7.6  x 1010  cm-3 

n (He)  = 3.4  x |07  cm"3  . 

For  hydrogen  we  also  follow  Jacchia  and  take  the  concentration  at  500km  to  vary 
with  according  to  the  relation 

log10  n(H;  500)  = 73.13  - 39.40  log,0  + 5.5  (log10  T J2  . 

The  thermal  diffusion  factor  a is  taken  as  -0.4  for  helium  and  as  zero  for  the 
other  constituents. 
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SUBROUTINE  ORBINT 

Summary  - The  subroutine  starts  and  controls  the  integration  of 

an  elliptic  orbit.  During  the  integration,  which  may 
be  forwards  or  backwards  in  time,  it  initiates  inter- 
polations at  constant  time  intervals  and,  if  required, 
will  determine  the  apses  of  the  orbit.  The  time  of 
the  apse,  the  osculating  elements  and  radial  distance 
are  written  to  an  unformatted  disc  file. 


Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (February  1975). 

Subroutine  statement  - SUBROUTINE  ORBINT. 

Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CINTEG/, 

/CONST/,  /CSMOON/  and  /INTP/  are  used  as  follows: 


Input  arguments  in  /CINTEG/  - 


MJDOCH 

N(I) 

PP  (3) 


PDV ( 1 8) 
PDVEL ( I 8) 


Modified  Julian  day  number  of  epoch. 

The  number  of  equations  to  be  integrated. 

Satellite's  geocentric  position  components, 

£ = (x,y,z)  , initially  at  epoch  (r^)  and  subsequently 
at  the  latest  integration  step  (MJDT + TIMET)  (km). 

Partial  derivatives  of  position  at  the  latest  integra- 
tion step  a_r/8(r0,vQ)  . 

Derivatives  with  respect  to  the  independent  variable 
s of  the  partial  derivatives  of  position 

_d_ 
ds 


ar 


^O’V 


V (3) 


A(  3) 


Satellite's  geocentric  velocity  components, 
v = (x,y,z)  , initially  at  epoch  (v^)  and  subsequently 
at  the  latest  integration  step  (km7s) . 

Satellite's  geocentric  acceleration  components, 

£ * (x.y.z)  at  the  latest  integration  step  (km/s  ). 

Partial  derivatives  of  acceleration  at  the  latest 
integration  step  ar/a(rQ,v0)  . 


PDACCT (18) 
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Input  arguments  in  /CONST/  - 

3 2 

EMU  Earth’s  gravitational  constant  (km  /s  ). 

Input  arguments  in  /CSM00N/  - 

MJDT  Modified  Julian  day  number  of  the  current  time. 

TIMET  The  current  time,  in  fractions  of  a day,  relative  to 

MJDT. 

Input  arguments  in  /INTP/  - 

EP  The  time  of  epoch  measured  in  fractions  of  a day 

relative  to  MJDOCH. 

TINT  The  time  at  which  the  next  interpolation  is  required. 

Initially  this  is  in  hours  from  epoch  and  subsequently 
in  days  from  epoch.  If  the  integration  is  backwards 
in  time,  TINT  must  be  negative. 

TIOR  Duration  of  the  integration  (hours) . 

Output  arguments  in  /CINTEG/  - 

XVEL,  YVEL,  ZVEL  The  initial  values  of  dx/ds,  dy/ds,  dz/ds  . 

T The  time,  in  seconds,  of  epoch  relative  to  MJDOCH 

(i.e.  T = EP  « 86400.0). 

Input  and  output  argument  in  /CINTEG/  - 

TVEL  The  value  of  dt/ds  at  the  latest  integration  step 

(s  *)  . 

Input  and  output  arguments  in  /INTP/  - 

DT  On  input,  the  interpolation  interval  in  minutes.  On 

output,  the  interpolation  interval  in  fractions  of  a 
day.  (If  the  integration  is  backwards  in  time,  DT 
must  be  negative.) 

NAP  Flag  controlling  the  detection  of  apses.  Initially 

NAP  = 1 if  apses  are  to  be  found  and  subsequently 
NAP  = N + 1 where  N is  the  number  of  apses  already 
found. 

Local  storage  used  - 69  real  array  variables,  17  real  variables  and  1 integer 

variable. 
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Subordinate  subprograms  - The  subroutines  COTOEL,  DEQRSPT,  INTPTA,  INTPTB,  MODAT, 

OUTPUT,  POAUXP,  SETPD,  SMPOS  and  TRINV  and  the 
functions  ANGLE , ANGL1 , INTFRC,  SOLVIN  and  SCPROD. 

Source  deck  - 60  cards  (ICL  code). 

Exp lanation  - The  subroutine  starts  and  controls  the  forward,  or 

backward,  integration  of  a satellite  orbit.  The  variables  TVEL  (dt/ds) , XVEL 
(dx/ds),  YVEL  (dy/ds)  and  ZVEL  (dz/ds)  are  set  in  terms  of  the  independent 
variable  s . If  a covariance  matrix  is  to  be  propagated,  the  subroutine  SETPD 
is  called  to  set  initial  values  of  the  array  PDVEL.  T is  set  to  the  time  of 
epoch  in  seconds,  and  TIOR  and  DT  are  converted  to  fractions  of  a day.  TINT,  the 
time  of  the  first  interpolation  relative  to  epoch  is  input  with  units  of  hours 
and  is  converted  into  a time  in  days  and  fractions  of  a day. 

The  integration  is  started  by  calling  DEQRSPT  with  the  statement  CALL 
DEQRSPT  (IND,  PDAUXP) . On  the  first  call  IND  = 1 and  on  the  second  and  subse- 
quent calls  IND  = 2 . 

After  each  integration  step,  the  time  elapsed  from  epoch  is  calculated  in 

days.  If  a table  of  apses  is  required,  the  length  of  the  radius  vector,  r.  , 

and  its  rate  of  change  r.  . are  determined.  If  r.r.  . < 0 there  is  an  apse 

l ’ l l-l  r 

between  the  last  two  integration  steps  and  SOLVIN  is  used  to  find  its  time  and 
radial  distance.  INTPTB  is  called  to  find  the  position  and  velocity  components 
at  the  apse  and  these  are  converted  to  osculating  elements  using  COTOEL.  This 
information  is  written  unformatted  to  disc  channel  4. 

A check  is  made  to  see  if  the  next  required  interpolation  falls  between 
the  last  two  integration  steps.  If  so,  INPTA  is  called  to  carry  out  this  inter- 
polation and  any  others  which  fall  between  the  two  steps.  When  this  is  completed, 
INTPTB  sets  TINT  to  the  time  of  the  next  required  interpolation  and  returns 
control . 

Before  the  next  integration  step,  the  geocentric  position  and  velocity 
components,  and  if  the  matrix  is  to  be  propagated,  the  partial  derivatives  of 
position,  velocity  and  acceleration  are  stored.  The  integration  is  continued 
until  the  time  of  the  latest  step,  measured  from  epoch,  exceeds  TIOR. 
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SUBROUTINE  OUTPUT 


- The  subroutine,  given  a set  of  earth  satellite  geocen- 
tric position  and  velocity  components  (POSVEL) , pro- 
vides any  combination  of  the  following  output  options: 

(i)  Write  POSVEL  to  line  printer. 

(ii)  Convert  POSVEL  to  osculating  elements  and  write 
to  disc. 

(iii)  Compute  range,  range  rate,  azimuth  and  elevation 
of  the  satellite  from  the  specified  ground 
stations  and  write  to  disc. 


Language 

Author 

Subroutine  statement 
Input  arguments 
PT(3) 


VT(3) 


MJD 

T 

P 

Use  of  COMMON 


- 1900  Fortran. 

- M.D.  Palmer  (February  1975). 

- SUBROUTINE  OUTPUT  (PT,  VT,  T,  MJD,  P) . 

Satellite's  geocentric  position  components  (x,y,z)  at 
epoch  MJD  + T (km) . 

Satellite's  geocentric  velocity  components  (x,y,z)  at 
epoch  MJD  + T (km) . 

Modified  Julian  day  number  of  epoch. 

Time,  in  fractions  of  a day  relative  to  MJD. 

2 2 2 I 

Radial  distance  at  epoch  (P  » (x  + y + z )2)  (km). 

- Certain  arguments  in  the  common  blocks  /STATION/, 
/CONST/  and  /INTP/  are  used  as  follows: 


Input  argument  in  /CONST/  - 

3 2 

EMU  Earth's  gravitational  constant  (km  /s  ). 

Input  argument  in  /INTP/  - 

NPV  If  NPV  > 0 , then  MJD,  T,  POSVEL,  radial  distance 

and  velocity  are  written  to  the  line  printer. 

Input  argument  in  /STATION/  - 


1 29 


IS 


Number  of  ground  stations. 


129 


Appendix  B 


53 


Input  and  output  arguments  in  /INTP/  - 

NEL  On  input,  if  NEL  > 0 , the  POSVEL  is  converted  to 

osculating  elements  and  written  to  disc.  If  this 
occurs  NEL  is  incremented  by  1. 


NAE 


Local  storage  used 
Source  deck 


On  input,  if  NAE  > 0 , the  range,  range  rate,  azimuth 
and  elevation  of  the  satellite  relative  to  the  speci- 
fied ground  stations  are  calculated  and  written  to 
disc.  If  this  occurs,  NAE  is  incremented  by  1. 

- 10  real  variables,  1 integer  variable. 

- 29  cards  (ICL  code). 


Subordinate  subprograms  - The  subroutines  COTOEL,  RAZEL  and  TRINV  and  the 

function  ANGLE. 


Explanation  - The  subroutine  is  provided  with  an  epoch  and  a corres- 

ponding POSVEL  and  radial  distance.  A number  of  output  options  are  available 
and  these  are  controlled  by  the  flags  NPV,  NAE  and  NEL. 

If  NPV  > 0 , the  satellite's  velocity  is  calculated  and  then  MJD,  T, 
x,  y,  z,  radial  distance,  x,  y,  z and  velocity  are  written  to  the  line  printer. 

If  NEL  > 0 , the  POSVEL  is  converted  to  osculating  elements,  using  sub- 
routine COTOEL.  MJD,  T , the  osculating  elements  and  radial  distance  are  then 
written  to  an  unformatted  disc  file  on  channel  3 and  NEL  incremented  by  1. 

If  NAE  > 0 , the  sidereal  angle  at  epoch  is  calculated  and  then  look  angles 
for  each  station  arc  found  using  subroutine  RAZEL.  For  each  ground  station, 

MJD,  T,  azimuth,  elevation,  range  and  range  rate  are  written  to  the  unformatted 
disc  file  on  channel  3 and  NAE  incremented  by  | . 
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SUBROUTINE  PDAUXP 

Summary  - The  subroutine  calculates  the  cartesian  components 

of  the  geocentric  accelerations  due  to  the  gravitational 
attractions  of  the  earth,  sun  and  moon,  atmospheric 
drag  and  the  precession  of  the  earth's  polar  axis.  It 
will  also,  if  required,  evaluate  the  partial  deriva- 
tives of  the  acceleration  components  with  respect  to  the 
initial  components  of  the  position  and  velocity  at 
epoch. 

Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (November  1974). 

Subroutine  statement  - Subroutine  PDAUXP. 

Use  of  COMMON  - Certain  quantities  in  the  conmon  blocks  /CINTEG/,  /CON/, 

/CONST/,  /CSM00N/  and  /PETURB/  are  used  as  follows: 


Input  arguments  in  /CINTEG/  - 


MJDOCH 

NEQ 

IAUX 
X,  Y,  Z 
PD( 18) 

XVEL,  YVEL,  ZVEL 
T 

PDVEL( 1 8) 


Modified  Julian  day  number  of  epoch. 

Number  of  equations  being  integrated  (22  if  partial 
derivatives  of  acceleration  are  required,  4 otherwise). 

Flag,  see  explanation. 

Cartesian  position  components  (x,y,z)  = r . 

Partial  derivatives  of  position  with  respect  to  the 
initial  position  (tq)  and  velocity  (Vq)  ; 3£./3(ro>vQ)  • 

The  latest  values  of  dx/ds,  dy/ds  and  dz/ds  . 

Time,  t , relative  to  MJDOCH  (seconds). 

Derivatives  with  respect  to  the  independent  variable 
s of  the  partial  derivatives  of  position 


8r 


¥(r0'V 


_d_ 

ds 
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Input  arguments  in  /CON/  - 


EJ 5 , EJ6,  EJ7,  EJ8,  EJ9 


C32 , S 32  1 
C41 , S41  > 
C43,  S43 J 


Coefficients  of  the  earth's  zonal  harmonics, 
J5’  J6 J9  • 


Fully  normalized  coefficients  of  certain  tesseral 


harmonics,  viz.  S32,  C^,  S^,  C43,  S43 


Input  arguments  in  /CONST/  - 


EMU 

ej2,  ej3,  ej4 


C22 , 

S22 

C31  , 

S3 1 

C33, 

S33 

C42 , 

S42 

C44 , 

S44 

ERAD 


E OMEGA 

SUNMU 

SELMU 


Earth's  gravitational  constant,  y 

Coefficients  of  the  earth's  zonal  harmonics,  J2,  J3, 
J4  ' 

Fully  normalized  coefficients  of  certain  tesseral 
harmonics,  viz.  ^"22*  ^22*  ^"31  * ^31  ' ^33'  ^33*  ^"42* 

S42 ' C44’  S44  - 

Mean  equatorial  radius  of  the  earth,  R . 

Mean  rotation  rate  of  the  earth's  polar  axis,  in  rad/s. 
Sun's  gravitational  constant,  yg  . 

Moon's  gravitational  constant,  y 

m 


Input  arguments  in  /CSMOON/  - 


XS,  YS,  ZS 


XM,  YM,  ZM 


Cartesian  components  of  the  sun's  position  (r  ) at 

s 

the  current  time,  computed  if  IAUX  < 0 . 

Cartesian  components  of  the  moon's  position  (r  ) at 

m 

the  current  time,  computed  if  IAUX  < 0 . 


Input  arguments  in  /PETURB/  - 


SMEAN 


SOL BAR 


RATE 

ARMA 


Mean  value  over  three  days  of  the  solar  10.7cm  radia- 

-22.  -2  -1 

cion  flux  in  units  of  10  “17  m Hz 

Mean  value  over  three  solar  rotations  of  the  10.7cm 

...  . -22  -2  -1 

radiation  flux  in  units  of  10  W m Hz  . 

Angular  velocity  of  the  atmosphere,  w in  rad/s. 

a 

The  product  AM  * HALFCD  x 1000  where  AM  is  the  satel- 
lite's area-to-mass  ratio  and  HALFCD  is  the  quantity 
$CD  where  CD  is  the  satellite's  drag  coefficient. 
If  ARMA  <=  0.0  , drag  terms  are  excluded. 
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TNITE 

LSP 

NSR 


Minimum  night-time  temperature  (K) . 

If  LSP  > 0 , luni-solar  perturbations  are  included. 

If  NSR  > 0 , solar  radiation  pressure  terms  are 
included  but  see  explanation. 


Output  arguments  in  /CINTEG/ 
XACC,  YACC,  ZACC 


2 2 2 2 2 2 

Values  of  d x/ds  , d y/ds  , d z/ds  at  the  current 


TVEL 

PDACC( 18) 


time. 

Value  of  dt/ds  at  the  current  time. 

Second  derivatives  with  respect  to  the  independent 
variable  s of  the  partial  derivatives  of  position 

j2 


ds 


31 


3(r0’V 


XVELT,  YVELT,  ZVELT  Cartesian  components  of  velocity,  x,  y,  z,  at  the 

current  time. 


XACCT,  YACCT,  ZACCT  Cartesian  components  of  acceleration,  x,  y,  z,  at  the 

current  time. 


PDACCT ( 1 8) 


Partial  derivatives  of  acceleration  with  respect  to 
the  initial  position  and  velocity  at  epoch 


dt 


3r 


8(r0’V 


Output  arguments  in  /CSMOON/ 


Modified  Julian  day  number  of  the  current  time. 

The  current  time,  in  fractions  of  a day,  relative  to 
MJDT. 

155  cards  (ICL  code). 

51  real  variables  and  3 integer  variables. 

The  subroutines  MODAT,  SMPOS  and  TRINV  and  the  functions 
ANGLE,  ANGL1,  INFIND  and  SCPROD. 


MJDT 
TIMET 

Source  deck 
Local  storage  used 
Subordinate  subprogram 

Explanation  - When  the  subroutine  is  called  for  the  first  time,  or  f r 

a change  of  epoch,  IAUX  must  be  set  to  -1.  For  subsequent  entries,  when  the 
time  has  changed  since  the  previous  entry,  IAUX  should  be  set  to  zero.  If  the 
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time  has  not  changed  IAUX  may  be  set  to  ! to  prevent  recomputation  of  certain 
terms.  (All  these  cases  are  allowed  for  by  the  calling  segment  in  POINT.) 

The  cartesian  components  of  the  geocentric  acceleration  acting  on  the 
satellite  are  found  by  adding  contributions  from  the  gravitational  attractions 
of  the  earth,  sun  and  moon,  atmospheric  drag  and  the  precession  of  the  earth's 
polar  axis.  The  subroutine  is  structured  so  as  to  permit  the  easy  insertion  of 
solar  radiation  pressure  at  a future  date. 

The  earth  harmonic  terms  are: 


u 

e 


I 

n ,ra 


I rn+2« 


, (m+ 1 ) 

n+ 1 


i) 


+ mW^e-jV^Cx  + j£) 


where  the  sidereal  angle  0 = 0^  + u^t  , 

W = (x  + jy)/rej0)  , 

„(ra) 


and 


is  the  mth  derivative  of  the  Legendre  polynomial,  P (z/r) 

n 


E0,0  “ 1 ’ 
En,0  = ' Jn 


2 = [2(2n  + 1 ) (n  - m)l/(n  + m)l]*(C  - jS  ) . 

n,ra  / J n,m  n,m 


Terms  are  included  up  to  the  ninth  zonal  harmonic  and  the  4,4  tesseral 
harmonic. 

The  accelerations  due  to  the  gravitational  attractions  of  the  sun  and 
moon  are  given  by 


b being  's'  for  the  sun  and  'm'  for  the  moon. 

Atmospheric  drag  is  included  if  ARMA  # 0 . A modified  version  of  the 
Jacchia  1965  model^  is  used  to  find  the  ambient  air  density.  Firstly  the  sub- 
routine computes  the  satellite's  height  and  latitude  (4>)  , the  declination  of  the 
sun  (6)  and  the  hour  angle  (H)  of  the  satellite  relative  to  the  sun.  The 
exospheric  temperature  (T^)  is  determined  using  the  equation. 
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T = T, 


N1TE 


1 + 0.280  + 0.28  [cos  - ej  [cos  jJ  j , 


where  0 = sin  L?  +.  j J. 


2.5 


and 


[h  - 0.78539816  + 0.20943951  sin  [l  + 0. 78539816] j . 


The  atmospheric  temperature  gradient  at  the  reference  altitude  (120km)  i« 
given  by: 


Tgrad0  = (Too  ~ 355X0.029  exp  [-  x2/2]) 

T^  - 800 

where  x = — - 

750  + 1.722  x |0_4(T  - 800) 2 

00 

Subroutine  MODAT  is  called  to  determine  the  atmospheric  density  p 
Ihe  force  acting  on  the  satellite  is  given  by  jp|v|2SCD  where  CD  is  the  drag 
coefficient,  S is  the  effective  cross-sectional  area  perpendicular  to  the  air 
flow  and.  V is  the  velocity  of  the  satellite  relative  to  the  ambient  air.  V is 
given  by 

V = v_  - a)  x £ 


where  _r  and  v are  the  position  and  velocity  vectors  of  the  satellite. 

If  the  cartesian  components  of  V are  V*,  Vy,  then  the  contributions 
to  the  acceleration  components  are: 

SCD 

x = - p 2T  lVlVx  * etc* 

where  M is  the  mass  of  the  satellite. 

SCL 

The  variable  ARMA  is  the  quantity  x 1000  , the  1000  being  a unit 

conversion  factor. 

The  contribution  to  the  acceleration  due  to  the  precession  term  is 


p(xz  - zx) 
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If  required,  the  subroutine  also  evaluate-?  the  partial  derivatives  of  the 
accelerations  with  respect  to  the  initial  position  (rQ)  and  velocity  (vQ)  at 
epoch,  i.e. 
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SUBROUTINE  POSVEL 

S umma ry  - The  subroutine  computes  the  geocentric  cartesian  com- 

ponents of  position  and  velocity  (POSVEL)  of  an  earth 
satellite  given  the  date,  time  and  a set  of  'PROP-type' 
orbital  elements. 


Language 

Author 

Subroutine  statement 


Input  arguments 
ELEMT(6,6) 


MJD 


TIME 


MJDOCH 


Output  arguments 
X,  Y , Z 

XDOT,  YDOT , ZDOT 
Use  of  COMMON 


- 1900  Fortran. 

- M.D.  Palmer  (February  1975). 

- SUBROUTINE  POSVEL  (ELEMT,  X,  Y,  Z,  XDOT,  YDOT,  ZDOT, 
MJD,  TIME,  MJDOCH). 


Array  of  coefficients  used  to  define  mean  orbital 
elements  (see  explanation). 

Modified  Julian  day  number  of  the  given  date. 

Time,  as  a fraction  of  a day  after  MJD,  at  which  the 
POSVEL  is  required. 

Modified  Julian  day  number  for  'PROP-type'  orbital 
elements. 

Note:  Normally  MJD  = MJDOCH  and  TIME  = 0.0  . 


POSVEL. 

- Certain  arguments  in  the  conmon  block  /CONST/  are  used 
as  follows: 


Input  arguments  in  /CONST/  - 


EMU 

EJ2 

ERAD 

Source  deck 
Local  storage  used 


Earth's  gravitational  constant. 

Earth's  second  zonal  harmonic,  J£  . 

Mean  equatorial  radius  of  the  earth,  R . 

- 115  cards  including  18  comment  cards  (ICL  code). 

- 4 integer  variables,  37  real  variables  and  60  real 
array  elements. 
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Subordinate  subprograms  - Subroutine  SHPAUX. 

Explanation  - The  subroutine  is  essentially  the  same  as  subroutine 

SATXYZ  used  in  the  program  PROP  ^ . It  will  normally  be  used  to  produce  a 
POSVEL  from  ' PROP- type' elements  at  the  epoch  corresponding  to  the  elements, 
i.e.  MJD  = MJDOCH  and  TIME  = 0.0  . The  standard  elements,  eccentricity  e , 
inclination,  i , right  ascension,  0 , argument  of  perigee,  u>  , and  mean 
anomaly,  M , are  assumed  to  be  polynomials  in  time  of  the  form 

2 3 4 5 

e = e.  + e,t  + e.t  + e_t  + e.t  + ect  . 

0 12  3 4 5 

The  coefficients  are  stored  in  the  array  ELEMT(6,6).  The  first  five  rows 
contain  six  coefficients  for  e,  i,  0,  u and  M and  the  first  five  elements  of 
the  sixth  row  are  set  to  the  values  of  the  corresponding  coefficients  in  the 
fourth-degree  polynomial  for  mean  motion. 

The  base  elements  (osculating  elements  with  short  and  long  periodic 
perturbations  removed)  are  computed  from  the  appropriate  polynomials,  whose 
degrees  are  specified  by  the  elements  of  the  array  N0MIAL(6),  set  in  a DATA 

statement.  For  the  standard  case  mentioned  above,  i.e.  t = 0.0  , the  elements 

of  NOMIAL  are  set  to  zero.  Long-periodic  perturbations  for  this  case  are  zero 
and  so  only  the  short-periodic  terms  need  be  incorporated  in  the  computation 
which  is  identical  to  that  in  PROP. 

- 2 1 

Firstly  a (semi-major  axis),  q (=  (1  - e )2),  E (eccentric  anomaly), 
2 _ 

v (true  anomaly),  p (aq  ) and  u (v  + 55)  are  computed.  Subroutine  SHPAUX 
then  computes  the  short  periodic  perturbations.  The  formulae  of  section  B.8  of 
Ref. 17  are  evaluated  to  obtain  the  POSVEL. 
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SUBROUTINE  PRINT 

Summary  - The  subroutine  reads  the  unformatted  disc  output  pro- 

duced by  subroutines  OUTPUT  and  ORBINT  and  writes  it 
to  the  line  printer. 


Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (January  1975). 

Subroutine  statement  - SUBROUTINE  PRINT  (IS,  NEL,  NAP,  NAE , K,  ST AN AM ) . 


Input  arguments 
IS 


NEL 


NAP 


NAE 


K 


STANAM (4 , 5) 


Use  of  COMMON 


Local  storage  used 


The  number  of  ground  stations  for  which  tables  of  look 
angles  have  been  produced. 

If  NEL  > 0,  NEL  = N + 1 where  N is  the  number  of 

lines  of  osculating  elements  stored  on  disc. 

If  NAP  > 0,  NAP  = N + I where  N is  the  number  of 

lines  in  the  table  of  apses  stored  on  disc. 

If  NAE  > 0,  NAE  = N + 1 where  N is  the  number  of 

lines  in  each  table  of  look  angles  stored  on  disc. 

The  number  of  different  sets  of  output  held  on  disc 
channel  3,  0 < K < 5 . 

Array,  holding  as  text,  the  names  of  the  ground 
stations  for  which  look  angles  have  been  determined. 

- None. 

- 10  real  variables  and  4 integer  variables. 


Subordinate  subprograms  - None. 


Source  deck 


- 50  cards  (ICL  code). 


Explanation  - The  subroutines  OUTPUT  and  ORBINT  store  unformatted 

information  on  disc  channels  3 and  4 respectively.  PRINT  recovers  this  informa- 
tion and  writes  it  to  the  line  printer  on  channel  2. 


The  information  produced  by  OUTPUT  falls  into  two  categories,  namely 
(a)  tables  of  time,  osculating  elements  and  radial  distance,  and  (b)  tables  of 
time,  azimuth,  elevation,  range  and  range  rate  for  up  to  four  ground  stations. 
The  data  is  stored  on  disc  channel  3 consecutively  and  in  chronological  order. 
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PRINT  recovers  and  lists  the  information  from  each  table  separately,  one 
line  at  a time.  The  table  of  osculating  elements  is  recovered  first  followed  by 
the  tables  of  look  angles.  For  the  latter,  only  those  lines  in  which  the  eleva- 
tion >0.0  are  actually  printed. 

ORBINT  produces  a table  of  apses  which  is  written  one  line  at  a time  to 
disc  channel  A.  PRINT  recovers  this  information  and  writes  it  directly  to  the 
line  printer.  If  no  apses  have  been  detected,  an  appropriate  comment  is  printed. 

A title  is  printed  before  each  table. 
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SUBROUTINE  RAZEL 

Summary  - The  subroutine  calculates  the  azimuth,  elevation,  range 

and  range  rate  of  an  earth  satellite  for  a given  ground 
s tation. 

Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (March  1975). 

Subroutine  statement  - SUBROUTINE  RAZEL  (SD,  I,  XX,  YY,  ZZ,  XD,  YD,  ZD,  AZ , 

EL,  R,  RR). 


Sidereal  angle  at  epoch  (rad). 

Station  number  in  the  range  1 c I < 4 . 

Satellite's  geocentric  position  components  (x,y,z)  (km). 

Satellite's  geocentric  velocity  components  (x,y,z) 

(km/s) . 


Azimuth  (rad). 

Elevation  (rad). 

Range  (km) . 

Range  rate  (km/s). 

- Certain  arguments  in  the  common  blocks  /STATION/  and 
/CONST/  are  used  as  follows: 

Input  arguments  in  /STATION/  - 


SL(4) 

Sin  , where  p.  is 

(i  = 1,4). 

the 

geodetic  station  latitude 

CL  (4) 

Cos  p.  . 

A1  (4) 

2 

R cos  4. 
e ri 

♦ h. 

l 

cos  4>.  where  R is 

i P 

/ 2 2 2.2 
R cos  4.  + R sin  4 . 
\ e l p l 

T 

the  earth's  polar  radius, 

R 

e 

is  the  earth's  mean 

equatorial  radius  and 

h. 

l 

is 

the  station's  height 

above  the  geoid  (i  “ 1 

,4). 

Input  arguments 
SD 
I 


Output  arguments 
AZ 
EL 
R 
RR 

Use  of  COMMON 


129 


Appendix  B 


65 


R sin  <}> . 

A2(4>  — ? — TT  + hi  sin  ^ = 1.4). 

^Re  cos  <t>^  + R^  sin  <fijJ 

S LONG (4)  The  geocentric  station  longitude,  A^  , rad  (i  = 1,4). 

Input  argument  in  /CONST/  - 

EOMEGA  Mean  rotation  rate  of  the  earth,  (rad/s). 

Subordinate  subprograms  - None. 

Source  deck  - 19  cards  (ICL  code). 

Local  storage  used  - 5 real  variables. 

Explanation  - POINT  will  produce,  as  an  output  option,  an  ephemeris 

of  look  angles  (azimuth,  elevation,  range  and  range  rate)  for  each  of  four 
ground  stations.  Each  call  to  RAZEL  gives  a set  of  look  angles  for  one  station 
for  a given  point  in  time,  the  station  being  identified  by  the  value  of  I 
(1  •<  I < 4). 

Consider  a cross-section  of  the  earth  containing  the  north  and  south  poles 
and  the  station.  Assume  an  axis  system  in  this  plane  so  that  0^  lies  along 
the  semi-major  axis  of  the  ellipse,  0 along  the  minor  axis  and  0 completes 
the  right-handed  set.  The  coordinates  of  a ground  station  are  given  by 

2 

R cos  <p 

X = — 5 9— r + h cos  <p 

^Re  cos  <t>  + Rp  sin  <f>J 


y = 0.0  • 

s 


z 

s 


R cos  (f 




(R 


2 2 , _2  2 . \ J 

cos  4>  + R sm  <J> ) ‘ 
n / 


♦ h sin  tf> 


where  <4  is  the  station’s  geodetic  latitude  and  h is  its  height  above  the 
spheroid.  The  satellite's  coordinates  (x,y,z)  in  PROP  axes  may  be  transformed 
into  the  above  axis  sytem  by  a rotation  about  the  z-axis  through  an  angle  of 
(A  ♦ 9),  where  A is  the  station's  geocentric  longitude  and  9 is  the  sidereal 
angle.  Thus  the  satellite's  coordinates  become 
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x'  = x cos  (A  + 0)  + y sin  (A  + 0) 
y'  = y cos  (A  + 0)  - x sin  (A  + 0) 


and  the  range  (r)  of  the  satellite  is  given  by 


■ [“ 


' - xg)2  + (y*  - yg)2  + (z’  - zg)2 


The  satellite's  transformed  velocity  components  are 


x'  = x cos  (A  + 0)  + y sin  (A  + 0)  + u^y 


y'  = y cos  (A  + 0)  - x sin  (A  + 0)  - ugZ 


where  u>  is  tne  earth's  rotation  rate.  The  satellite's  range  rate  (r)  is 
given  by 
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r = (x'x'  + y'y ' + z'z')/r 

The  azimuth  (A)  and  elevation  (E)  of  the  satellite  from  the  station  are 
found  in  a 'local'  axis  sytem  with  the  origin  at  the  station,  the  x-axis  along 
the  local  vertical,  the  z-axis  along  the  local  north  and  the  y-axi9  in  the 
direction  of  local  east.  The  coordinates  c^e  satellite  in  this 

system  are: 

= (x'  - x ) cos  $ ♦ (z'  - z ) sin  <t> 

s s 

= y’ 

= (z'  - z ) cos  $ - (x'  - x ) sin  $ 

s s 

A - ATAN2(y^,z^) 

E - ATAN2[xe,(y2  ♦ *J)*]  . 

A rigorous  derivation  of  this  theory  is  given  in  Ref. 18. 


Thus 

and 
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SUBROUTINE  READEL 

- The  subroutine  reads  RAE  orbital  elements  from  a set  of 
five  punched  cards  and  stores  them  in  a standard  array. 

- 1900  Fortran. 

- M.D.  Palmer  (March  1975). 

- SUBROUTINE  READEL  (ELEMT,  MJDOCH) . 


Array  of  coefficients  used  to  define  orbital  elements 
(see  explanation). 

MJD  of  the  midnight  epoch  at  which  the  elements  are 
defined. 

- None. 

- 15  cards  (ICL  code). 

- 6 real  array  elements,  3 real  variables  and  3 integer 
variables. 

Subordinate  subprograms  - None. 

Explanation  - The  subroutine  is  essentially  the  same  as  the  sub- 
routine ELREAD  used  in  the  program  PROP  it  is  assumed  that  the  five  basic 

orbital  elements  of  an  earth  satellite,  e,  i,  12 , w,  M may  be  represented  by 
polynomials  in  time  of  degree  five,  e.g. 

2 3 4 5 

e = e^  + ejt  + e2C  + e^t  + e^t  + e^t 

The  six  coefficients  for  each  polynomial  are  read  from  one  card  and  stored  in  one 
row  of  the  array  ELEMT.  The  first  five  elements  of  the  sixth  row  are  set  to  the 
values  of  the  corresponding  coefficients  in  the  fourth  degree  polynomial  in  mean 
motion  as  follows: 

ELEMT (6 , j ) = j ELEMT (5 , j + 1 ) for  j = 1,5  . 

The  coefficients,  as  read,  have  units  of  degrees  and  days  and  these  are 
changed  to  radians  and  seconds  before  storage. 

The  five  data  cards  must  contain  MJDOCH  in  cols. 1 to  5 and  a serial  number 
( 1 , . . . , 5)  in  coi.7.  A check  is  made  that  the  cards  are  in  the  correct  sequence 
and,  if  not,  a STOP99  instruction  is  obeyed. 


Summary 

Language 

Author 

Subroutine  statement 
Output  arguments 
EI.EMT  (6,6) 

MJDOCH 

Use  of  COMMON 
Source  deck 
Local  storage  used 


d 
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SUBROUTINE  REIP 

Summary  - The  subroutine  reads  the  input  data  for  POINT  and, 

where  necessary,  converts  it  to  the  form  required  by 
the  main  program. 

Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (January  1975). 

Subroutine  statement  - SUBROUTINE  REIP. 

Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CINTEG/, 

/CONST/,  /INTP/, /PETURB/  and  /STATION/  are  used  as 
follows: 


Input  arguments  in  /CONST/  - 

3 

EMU  Earth's  gravitational  constant,  y (km  /s) . 

ERAD  Mean  equatorial  radius  of  the  earth,  R (km). 

EOMEGA  Mean  rotation  rate  of  the  earth  w in  rad/s. 

e 

Input  arguments  in  /PETURB/  - 

HALFCD  The  quantity  $C^  where  C^  is  the  satellite's  drag 

coefficient. 


Output  arguments  in  /CINTEG/  - 


MJD 

NM 

X,  Y,  Z 


XD,  YD,  ZD 


Modified  Julian  day  number  of  epoch. 

Number  of  equations  to  be  integrated. 

Cartesian  components  of  the  satellite's  geocentric 
position  (x,y,z)  at  epoch  (km). 

Cartesian  components  of  the  satellite's  geocentric 
velocity  (x,y,z)  at  epoch  (km/s). 


Output  arguments  in  /INTP/  - 

EP  Time  of  epoch,  in  fractions  of  a day  relative  to  MJDOCH. 

DT  Interpolation  inteval,  min. 

TINT  The  time  from  epoch  at  which  the  first  interpolation  is 

required,  hours. 


TIOR 


Integration  period,  hours. 
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NPV 

NEL 

NAP 

NAE 

COV (6,6) 
NOTPT 


NPV  =1  if  an  ephemeris  of  geocentric  position  and 
velocity  components  (POSVEL)  is  required. 

NEL  =1  if  an  ephemeris  of  osculating  elements  is 
required. 

NAP  =1  if  a table  of  apses  is  required. 

NAE  = 1 if  tables  of  look  angles  are  required. 

Covariance  matrix  corresponding  to  the  POSVEL. 

NOTPT  =1  if  an  ephemeris  of  POSVEL  and,  if  specified 
the  propagated  covariance  matrix  is  to  be  written  to 
a permanent  disc  file. 


Output  arguments  in  /PETURB/  - 

ARMA  The  product  AM  * HAT.FCD  * 1000.0,  where  AM  is  the 

2 

area-to-mass  ratio  of  the  satellite  in  m /kg.  If 
ARMA  # 0.0,  drag  terms  are  included  in  the  integration 

TNITE  Nightime  minimum  exospheric  temperatures  (K.) . 

SMEAN  Mean  value  over  three  days  of  the  solar  10. 7cm 

. . -22  -2  -1 
radiation  flux  in  units  of  10  W m Hz 


SOLBAR 


RATE 

LSP 


NSR 


Mean  value  over  three  solar  rotations  of  the  10. 7cm 

. . . . -22  -2  -1 
radiation  flux  in  units  of  10  w m Hz 

Angular  velocity  of  the  atmosphere,  rad/s. 

LSP  = 1 if  luni-solar  perturbations  are  included  in 
the  integration. 

A spare  variable,  to  be  used  as  a flag  should  solar 
radiation  pressure  terms  be  incorporated  in  the 
integration. 


Output  arguments  in  /STATION/  - 

K The  number  of  ground  stations  for  which  tables  of 

look  angles  are  required  (K  ^ 4). 

SLT(4)  Sin  , where  <p . is  the  geodetic  station  latitude 

(i  = 1 , . . . ,4) . 

Cos  p . ( i = 1 

l 


CLT (4) 


» • • • » 


4). 
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Al  (4) 


A2(4) 

SLO  C4) 
STANAM(4 , 5) 


2 

R co  s d> . 
e i 


(R 


cos'-*.  + R“  s in- . J " 
e l p 1/ 


+ h.  cos  <p . where  R and 
li  p 


are  the  earth's  polar  and  mean  equatorial  radii 


and  Ik  is  the  station's  height  above  the  geoid 
(i  = 1 , . . . ,4) . 


R sin 
P i 


/ 2 


(R  cos  <p. 

\ e l p 


+ R2  sin2^.)^ 


+ h.  sin 


The  geocentric  station  longitude  X^  , rad  (i  = 1,, 


..4), 


Array  holding  the  ground  station  names  as  text. 


Input  and  output  arguments  in  /CINTEG/  - 

H On  input,  the  integration  step  length  as  set  in  BLOCK 

DATA.  On  output  H/Cl  (Cl  =£  0)  where  Cl  is  an  input 
data  item. 

Subordinate  subprograms  - The  subroutines  ELTOCO,  IPTOCO,  MATMUL,  POSVEL,  RLEASE , 

ROTATE,  SDC2EL,  SETPD , SHPAUX  and  TRAMAT  and  the 
functions  ANGLE  and  EAFKEP . 


Source  deck  - 130  cards  (ICL  code). 

Local  storage  used  - 113  real  array  elements,  18  real  variables  and  2 integer 

variables. 


Explanation  - The  subroutine  determines  the  position  and  velocity  of 

an  earth  satellite  at  a given  epoch  in  terms  of  its  geocentric  cartesian  compo- 
nents in  the  PROP  axes  system.  The  components  may  be  read  directly  or  input  in 
an  alternative  form  and  converted.  The  four  alternative  forms  are: 


(i)  a standard  set  of  launch  vehicle  injection  conditions  (speed,  climb  angle, 
azimuth,  radius,  latitude  and  longitude); 

(ii)  a set  of  osculating  elements  (semi-major  axis,  eccentricity,  inclination, 
right  ascension  of  the  ascending  node,  argument  of  perigee  and  mean 
anoma ly) ; 

(iii)  a set  of  USAF  Space  Defense  Center  (SDC)  two  line-elements;  and 

(iv)  a set  of  RAE  five-line  mean  elements. 
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A choice  of  the  perturbations  to  be  included  in  the  integration  is  provided 
and  there  are  a number  of  output  formats. 

One  option  is  to  propagate  the  covariance  matrix  associated  with  the 
position  and  velocity  components.  The  initial  covariance  matrix  may  be  provided 
either 

(i)  in  terms  of  the  position  and  velocity  components  (km  and  km/s) 

or  (ii)  in  terms  of  injection  conditions  with  units  of  ft/s,  degrees  and 
feet.  In  this  case,  the  orbit  must  also  be  supplied  as  a set  of 
injection  conditions. 

In  the  latter  case,  the  subroutine  carries  out  a transformation  of  the 


cov(Xq)  = A cov(Yq)  AT 


where  eov (X^j  is  the  covariance  matrix  in  terms  of  position  and  velocity 
components , 

cov  (Yq)  is  the  injection  covariance  matrix 
and  A is  the  matrix  of  pa'cial  derivatives  given  by 
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where  x,  y,  z are  the  position  components, 
x,  y,  z the  velocity  components, 

v the  speed  at  injection, 

0 the  climb  angle, 

<1/  the  azimuth, 


i 
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and 


0 is 


where 


Row  1 


Row  2 


Row  3 


Row  4 


X’ 


the  radial  distance, 

the  geocentric  latitude 

the  longitude  in  inertial  space. 


X'  is 

defined  by  X'  = 0 + X 

where 

X is  the  geocentric  long: 

the  modified  sidereal  angle  given  by 

0 = 

100.075542  + 

360.9856 1 2 288 ( t - 33282.0) 

t = MJD  + EP  . 

It  can 

be  shown 

that: 

3 x 
3v 

9x  9x 

90  ~ 9^ 

0.0  , 

9x  , . , 

= cos  4 cos  X 
9r 

3x 

34> 

- r sin  <(>  cos 

X’  , 

9x  , « . f 

-rrr  = - r cos  41  sin  X 

a A 
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3v 

9y  _ 9y 
90  = dip 
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= cos  ip  sin  X' 
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iz  = 

34) 

- r sin  <p  sin 
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= r cos  p cos  X' 
o A 
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9 z _ 9 z 

90  = 9* 
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9z 

— = sin  * 
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■ sin  0 
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Row  5 


t 


Row  6 


3x 

Jp 

3x 

JF~ 


ii  = 

3v 

ii  = 

39 

3y  , 

3p 

2i  = 

3 r 

ii  = 

3 d> 

_ 

slT 


v[-  sin  0 sin  <J>  cos  A'  - cos  0 cos  p cos  p cos  A 


v [-  sin  0 cos  4>  sin  A'  + cos  0 [cos  p sin  p sin  A'  - sin  p cos  A']] 


sin  0 cos  (J>  sin  A'  + cos  0 [sin  d>  cos  A'  - cos  p sin  <p  sin  A'} 


v [cos  0 cos  <}>  sin  A'  - sin  0 [sin  ip  cos  A'  - cos  p sin  p sin  A']] 


v cos  0 [cos  ip  cos  A'  + sin  4>  s in  ip  sin  A'] 
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Summary 

Language 

Author 

Subroutine  statement 
Input  arguments 
ANGLE 

ANCVEL 


SUBROUTINE  ROTATE 

- The  subroutine  performs  rotations,  about  the  z-axis 
by  a given  angle  at  a given  angular  velocity. 

- ASA  Fortran  (Standard  Fortran  4). 

- A.W.  Odell  (August  1970). 

- SUBROUTINE  ROTATE  (X,  Y,  XDOT,  YDOT,  ANGLE,  ANGVEL) . 


Angle  B , through  which  coordinates  are  to  be 
rotated  - anticlockwise. 

Angular  velocity,  <o  , to  be  applied  to  velocity 
coordinates  - anticlockwise. 


Input  and  output  arguments  - 


X 

Y 

XDOT 

YDOT 

Use  of  COMMON 
Source  deck 


x coordinate, 
y coordinate, 
x coordinate, 
y coordinate. 

- None. 

- 11  cards  (ICL  code). 

- 3 real  variables. 


Local  storage  used 
Subordinate  subprograms  - None. 

Explanation  - The  rotation  is  equivalent  to  the  following 


(x,y)  -*•  (x  cos  9 - y sin  9 , y cos  0 + x sin  9) 


(x,y)  •+ 


followed  by 


(x  cos  9 - y sin  9 - uy  , y cos  9 + x sin  9 + wx) 
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Summary 


FUNCTION  SCPROD 

- The  function  gives  the  scalar  (inner)  product  of  two 
arrays . 


Language 

Author 

Function  statement 
Input  arguments 
A(l),  B ( 1 ) 


- PLAN,  for  use  with  1900  Fortran. 

- A. W.  Odell  (March  1973)  . 

- FUNCTION  SCPROD  (A,  B,  NA,  NB,  N) . 


Locations  of  the  first  elements  to  be  multiplied  (must 
be  array  elements). 


NA,  NB 
N 


Increments  in  the  arrays.  A,  B of  elements  to  be 
multiplied. 

Number  of  elements  to  be  multiplied. 


Output  function 
SCPROD 

Use  of  COMMON 
Source  deck 
Local  storage  used 


The  scalar  product 


N 

^ A ( 1 + (i  - 1 )NA) B( 1 + (i  - 1 )NB) . 
i = l 


- None. 

- 25  cards,  including  1 conment  card  (ICL  code). 

- 2 words  for  variables,  plus  22  words  for  program. 


Subordinate  subprograms  - None. 

Explanation  - The  arrays  in  the  calling  routine  may  have  any  dimen- 

sions although  they  are  treated  as  one-dimensional  arrays  in  the  function.  For 
simplicity  A and  B are  dimensioned  '1'.  This  may  cause  the  function  to  fail 
in  TRACE2,  so  it  should  not  be  compiled  in  this  mode.  If  N ^ 0 , the  function 
will  return  zero. 
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SUBROUTINE  SDC2EL 

Summary 

Language 
Author 

Subroutine  statement 
Output  arguments 
MJDOCH 
ELEMT (6,6) 

Use  of  COMMON 

Input  arguments  in  /CONST/  - 

3 2 

EMU  Earth's  gravitational  constant,  y (km  /s  ). 

XJ2  Earth's  second  zonal  harmonic  coefficient,  . 

XJ3  Earth's  third  zonal  harmonic  coefficient,  J^  . 

ERAD  Earth’s  mean  equatorial  radius,  R (km). 

Subordinate  subprograms  - None. 

Local  storage  used  - 5 integer  variables  and  16  real  variables. 

Source  deck  - 71  cards  including  10  comment  cards  (ICL  code). 

Explanation  - This  subroutine  is  essentially  the  subroutine  SDCELS 

used  in  the  program  PROP  . A set  of  orbital  elements  consists  of  eccentricity, 
e , inclination,  i , right  ascension  of  ascending  node,  fl  , argument  of  perigee, 
u>  , mean  anomaly,  M , and  mean  motion,  n . With  the  exception  of  mean  motion, 
these  elements  may  be  assumed  to  be  polynomials  in  time,  i.e. 

2 

e - eQ  + ejt  + e2t  + ...  etc., 
where  e^ , e^,  etc.  are  the  'rate  elements'. 


- The  subroutine  reads  a set  of  Space  Defense  Center 
(SDC)  two-line  elements  and  converts  them  to  an 
equivalent  set  of  'PROP'  type  elements. 

- 1900  Fortran. 

- M.D.  Palmer  (December  1974). 

- SUBROUTINE  SDC2EL  (MJDOCH,  ELEMT). 

Modified  Julian  day  number  of  epoch. 

Array  of  PROP-type  orbital  elements. 

- Certain  arguments  in  the  common  block  /CONST/  are  used 
as  follows: 
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SDC  elements  differ  from  PROP  elements  in  the  following  ways: 

(i)  PROP  elements  are  for  midnight  epochs,  whereas  SDC  elements  relate 
to  an  ascending  node. 

(ii)  PROP  elements  follow  the  long-periodic  motion  of  the  satellite,  but 
SDC  elements  have  been  'corrected'  for  long-periodic  perturbations 
associated  with  a nominal  value  of  the  earth's  third  zonal  harmonic. 

(iii)  The  origin  of  the  right  ascension  of  the  ascending  node  is  the  mean 
equinox  of  date  for  SDC  but  the  1950.0  equinox  for  PROP. 

(iv)  SDC  elements  i q , fig,  Wq  are  in  degrees  and  the  SDC  units  for  the 
M-polynomial  quantities  are  revolutions  and  days,  whereas  PROP 
elements  are  stored  in  units  of  radians  and  seconds. 

The  quantities  read  from  the  first  SDC  card  are  card  number,  satellite 
identification  number,  year  of  epoch,  day  of  year,  fraction  of  day,  M?  and  . 

It  is  read  in  the  format 

(II,  IX,  15,  2X,  A8,  IX,  12,  13,  F9.8,  IX,  F10.8,  IX,  E8.5)  . 

The  second  card  contains  card  number,  satellite,  iQ,  0g,  e^,  Uq  , Mg  and 
M|  and  is  read  in  the  format 

(II,  IX,  15,  2( IX,  F8.4),  IX,  F7.7,  2(1X,  F8.4),  IX,  FI  1 .8)  . 

The  semi-major  axis,  which  is  an  auxiliary  element,  is  computed  from 
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The  time  in  days,  x , by  which  the  new  epoch  is  later  than  the  epoch  of 
the  SDC  elements  is  found.  fig,  uig  and  Mg  are  replaced  by  Og  + OjT,  Uq  + to  ^ T 
and  Mq  + MjT  + , respectively,  and  Mj  and  M2  by 

Mj  + 2M2T  + and  M2  + 3M^t  . Coefficients  of  the  mean-motion  polynomial 

are  obtained  from  n.  = (j  + l)Mj+j  for  j = 0,  1 , 2 and  3. 

fig  is  decreased  by  3.506  * 1 0 where  D is  the  number  of  days  elapsed 
from  1950.0,  to  allow  for  equinox  conversion. 

The  SDC  long-periodic  perturbations  are  computed  from 

K, 


ll  |~3  + 5c~ 

lL  4 [l  + c 


eQ  cos  w0 


and 


yNL 


K3/2 


where  c 
and 


K3  = J3  sin  i0/J2a(t  - eQ)  . 


Mg  is  replaced  by  Mg  + L^,  eg  by  [C  + S ] 


where  C 
M0  + “0 


eg  cos  Wg  and 
remains  constant. 


S = eQ  sin  ^ ♦ ayfJL 


and  o)g  by 


tan  1 (S/C)  , 


is  adjusted  such  that 


The  derived  PROP  elements  are  stored  in  the  following  positions  of  the 
array  ELEMT:  eg,  i Wg,  Mg  and  ng  make  up  the  first  column;  Mg  , in  the 

in  the  sixth  row  is  followed  by 


fifth  row,  is  followed  by  Mj , M3  and  M3 


0 


n^  and  n^ 


All  other  elements  of  ELEMT  are  set  to  zero. 


J 
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SUBROUTINE  SETPD 


Summary 


Language 

Author 


- The  subroutine  sets  the 

the  partial  derivatives 

_d_ 
ds 


3r 


9(vVj 


initial,  non-zero  values  of 
3r/3 (rgi Vq)  and 


- 1900  Fortran. 


- M.D.  Palmer  (February  1973). 


Subroutine  statement  - SUBROUTINE  SETPD. 

Use  of  COMMON  - Certain  arguments  in  the  common  block  /CINTEG/  are  used 

as  follows: 


Input  argument 
TVEL 


The  initial  value  of  dt/ds  . 


Output  arguments 
PD  ( 1 8) 

PDVEL( 18) 

Source  deck 

Local  storage  used 

Subordinate  subprograms 


Array  holding  the  partial  derivatives 
3_r /3  (rQ»vQ}  • 

Array  holding  partial  derivatives 


of  position 


3r 


8 cards  (ICL  code) . 
1 integer  variable. 
None. 


Explanation  - Before  using  DEQRSPT  to  integrate  the  22  equations  >f 

motion,  the  initial  non-zero  values  of  two  partial  derivative  matrices  must  be 
set  as  follows: 


3r 


3(ro’vo} 


PD ( 1 ) ...  PD(6) 

PD(7)  ...  PD (12) 

PD ( 1 3)  ..  PD ( 1 8) 


3x 

3x 

3x 

3x 

3x 

3x 

3xq 

3yo 

3zo 

3x0 

3yo 

3i0 

= 

iz_ 

3xq 

iiL. 

3yo 

3^_ 

3zo 

h- 

3xq 

3yo 

i)L 

3io 

3 z 

3 z 

3 z 

3z 

8 z 

3 z 

_3x0 

3yo 

3zq 

3x0 

3*o 

3V 
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I 0 0 
0 1 0 
.0  0 1 


0 0 0 
0 0 0 
0 0 0. 


and 


'PDVEL(I)  ...  PDVEL(6)  " 

3£ 

PDVEL(7)  ...  PDVEL( 1 2) 

~di 

3(vV_ 

_ds_ 

PDVEL( 1 3)  ...  PDVEL( 1 8) 

3x 

3x 

3x 

3x 

3x 

3x 

3xq 

3yo 

3zo 

3x0 

3y0 

3zo 

ii_ 

ii_ 

ii- 

1L. 

iz_ 

3xq 

3yo 

3zo 

3xo 

3yo 

3zo 

3 z 

3z 

3 z 

3z 

3 z 

Si 

_3xo 

3yo 

3zq 

3xo 

3yo 

3zo. 

0 

0 

0 

Jt/ds 

0 

0 

M 

0 

0 

0 

0 

dt/ds 

0 

_0 

0 

0 

0 

0 

dt/ dsj 

where  s is  the  independent  variable  given  initially  by  dt/ds  = r^  p , 
r0  * ^x0,y0,Z0^  being  the  initial  radius  vector  to  the  satellite  and, 

Vq  = (xQ,yQ,Zg)  the  initial  velocity  of  the  satellite  and  p the  earth's  gravi- 
tational constant. 
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Summary 

Language 

Author 

Subroutine  statement 

Input  arguments 
ECC 
SI 
Cl 

OMEGA 

V 

SV 

CV 

U 

FMANOM 

F 

P 

POVERR 

Q 

EJ2 

ERAD 

EMU 

Output  arguments 
TSP(6) 

Use  of  common 


SUBROUTINE  SHPAUX 

- The  subroutine  calculates  the  short-periodic  perturba- 
tions, due  to  the  zonal  harmonic  , in  an  earth 

satellite's  orbit. 

- 1900  Fortran. 

- M.D.  Palmer  (January  1975). 

- SUBROUTINE  SHPAUX  (ECC,  SI,  Cl,  OMEGA,  V,  SV,  CV, 

U,  FMANOM,  F,  P,  POVERR,  Q,  TSP,  EJ2,  ERAD,  EMU). 


Eccentricity,  e . 

Sin  i , where  i is  the  inclin  cion. 

Cos  i . 

Argument  of  perigee,  ai  (rad). 

True  anomaly,  v (rad). 

Sin  v . 

Cos  v . 

Argument  of  latitude,  u (rad). 

Fractional  part  of  mean  anomaly,  (rad). 

2. 

Sin  l . 

Semi-l’atus  rectum,  p (km)  . 

p/r  , where  r is  the  satellite's  radial  distance  (km). 
(1  -eV  . 

Earth's  second  zonal  harmonic,  . 

Mean  equatorial  radius  of  the  earth  (km). 

3 2 

Earth's  gravitational  constant,  y (km  /s'”). 

Array  of  short-periodic  perturbations,  (see 
explanation) . 

- None. 
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Source  deck  - 30  cards  including  3 comment  cards  (ICL  code). 

Local  storage  used  - 14  real  variables. 

Subordinate  subprograms  - None. 

Explanation  - The  subroutine  is  essentially  the  same  as  the  sub- 

routine SHOPER  used  in  the  program  PROP  ^ , It  computes  the  short-periodic 

perturbations  required  by  the  subroutine  POSVEL.  Input  elements  are  assumed  to 

1 9 

be  the  mean  elements  of  Merson  , and  the  expressions  for  the  short-periodic 
perturbations  are: 

di  = iKC  sin  i cos  i , 
s ’ 

- K(v  - + e sin  v - JS)  sin  i cos  i , 

Kh  [v  - + e sin  v + 5 ge  sin  v (g  + cos  v)] 

+ jKf[i  sin  2u  + e sin  (2u  - v)]  , 

dpg/2p  = jK(hq  + fc)  , 
clrs/p  = lK[i  f cos  2u  - h ( 1 + g cos  v - qr/p)] 

and 

drg  = j(u/p)^K[h  sin  v (-  Jqe  + g(p/r)2)  - f(p/r)2  sin  2u] 
where  K = 1 iJ2(R/p)2, 

C = cos  2u  + e [cos  (2u  - v)  + j cos  (2u  + v)] , 

S = sin  2u  + e [sin  (2u  - v)  ♦ 5 sin  (2u  + v)], 
h = 1 - 1 i f 
and  g » e/( I + q)  . 

The  perturbations  are  stored,  in  the  above  order,  in  the  array  TSP. 


dO  sin  i = 
s 


du  + dO  cos  i = 
s s 
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Summary 

Language 

Author 

Subroutine  statement 
Dummy  arguments 
Use  of  COMMON 


SUBROUTINE  SMFOS 

- The  subroutine  computes  the  geocentric  cartesian 
coordinates  (km)  of  the  sun  and  moon  at  a given  time. 

- ASA  Fortran  (Standard  Fortran  4). 

- A. W.  Odell  (May  1973) . 

- SUBROUTINE  SMPOS. 

- None. 

- Certain  quantities  in  common  block  /CSMOON/  are  used 
as  follows: 


Input  arguments  in  /CSMOON/  - 

MJDT  Modified  Julian  day  number  of  the  required  time. 

TIMET  Time,  as  a fraction  of  a day,  since  0 hours  on  MJDT. 


The  cartesian  coordinates  (x  ,y  ,z  ) = r of  the  moon. 

m’-'m’  m m 


Output  arguments  in  /CSMOON/  - 

XS,  YS,  ZS  The  cartesian  coordinates  (x  ,y  ,z  ) = r of  the  sun. 

s s s s 

XM,  YM,  ZM 

Input  and  output  argument  in  /CSMOON/  - 
TABLE(43)  Working  space. 

Source  deck  - 27  cards  including  1 comment  card  (ICL  code). 

Local  storage  used  - 2 integer  variables,  1 logical  variable. 

Subordinate  subprograms  - The  subroutine  UTDA , and  the  functions  INFIND  and 

SCPROD. 


Explanation  - The  subroutine  obtains  the  sun-moon  coordinates  by 

interpolation  in  a table  of  daily  positions,  with  2nd  and  4th  differences,  using 
Everett's  interpolation  formula.  The  table  is  read  from  a disc  file  when 
required,  using  the  subroutine  UTD4.  It  is  stored  in  an  area  called  SUNMOONTABLE , 
the  location  of  which  is  found  using  the  function  INFIND.  The  area  is  contained 
in  a disc  file,  which  must  have  been  opened,  before  entry  to  SMPOS,  as  peripheral 
unit  ED 7. 


The  table  is  stored  as  follows:  firstly,  the  modified  Julian  day  numbers 

' the  first  and  last  sets  of  coordinates  (2  integers),  then  data  for  each 

2 2 4 4 

ught  as  follows  (sets  of  18  reals):  r,r,Ar,Ar,Ar,Ar 

s m’  s’  m’  s’  m 
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If  f„,  f denote  values  of  f at  t = 0 and  t = I , Everett's  inter- 

U . 20 
polation  formula  to  the  4th  differences “ gives: 

f(t)  = DQ  + th>  + (t  - 1)^D2  + (t  + 1)^D3  + (t  - 2)(^Da  + (t  + 2)d5))J)  , 

where  = 62nfQ/(2n)!,  D2n+]  = 62n(fj  - fQ)/(2n  ♦ 1)!. 

3 A 

The  data  was  originally  obtained  on  punched  cards  from  JPL  ’ and  has  been 
transformed,  before  storing  on  the  disc,  into  the  SA0/PR0P  ^ system  of  axes, 
using  the  subroutine  AX1950.  If  the  time  input  to  SMPOS  lies  outside  the  range 
of  data  stored  on  the  disc,  a STOP77  statement  will  be  obeyed. 
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FUNCTION  SOLVIN 


Summary 


Language 

Autnor 

Function  statement 
Input  arguments 
MODE 
H 

FO,  FI 
FDO , FD1 

F 

Output  function 


- Given  a function  and  its  derivative  at  two  points,  the 
function  solves  four  problems  involving  cubic 
interpolation 

- ASA  Fortran  (Standard  Fortran  4) . 

- A.W.  Odell  (January  1975). 

- FUNCTION  SOLVIN  (MODE,  M,  FO,  FI,  FDO,  FD1 , F)  . 

Number  specifying  problem  to  be  solved,  see  explanation. 

Difference,  h , between  the  arguments  of  the  function 
(i.e.  h = X]  - xQ) . 

Function  values  fQ  at  x = xQ  and  f ^ at  x = Xj  . 

Derivative  values  f^  , at  x = xQ  and  fj  at 
x = x,  . 

Required  function  value  or  argument  - see  explanation. 


SOLVIN 

Use  of  COMMON 
Source  deck 


Solution  to  problem  - see  explanation. 

- None. 

- 33  cards,  including  5 comment  cards  (ICL  code). 

- 9 real  variables. 


Local  storage  used 
Subordinate  subprograms  - None. 

Explanation  - A cubic  polynomial  P(x)  = f^  + f^x  + a7x~  + a^x2  is 

first  fitted  to  the  data  giving  a2  = (3B  - A)/h  and  a^  = (A  - 2B)/h2  where 
A = f'  - fQ  and  B = (f , - fQ)/h  - t'Q  . 


Then 

>> 

(1) 

if 

MODE  = 

1 

, Newton's  method  is  used 

to  find 

the 

value 

of  X 

for 

wh  i ch 

P(x) 

= F , 

(2) 

if 

MODE  = 

2 

, Newton's  method  is  used 

to  find 

the 

value 

of  X 

for 

which 

P'  (x 

) - F , 

(3) 

if 

MODE  = 

3 

, P (F)  is  evaluated,  and 

(4) 

if 

MODE  = 

4 

, P'(F)  is  evaluated. 
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I 


Sunmary 

Language 

Author 

Subroutine  statement 
Input  arguments 
A 
K 
L 


SUBROUTINE  TRAMAT 

- The  subroutine  performs  matrix  transposition. 

- 1900  Fortran. 

- M.D.  Palmer  (May  1973). 

- SUBROUTINE  TRAMAT  (A,  AT,  K,  L) . 

The  matrix  to  be  transposed. 

Number  of  rows  in  A. 

Number  of  columns  in  A. 


Output  arguments 
AT 

Use  of  COMMON 
Source  deck 


Local  storage  used 
Subordinate  subprograms  - None. 


The  transposed  matrix. 

- None. 

- 7 cards  (ICL  code). 

- 2 integer  variables. 


Explanation 


- The  subroutine  interchanges  the  rows  and  columns  of  a 


matrix  A to  form  the  transpose  A , i.e. 


AT(i,j)  = A(j,i) 


129 


Appendix  B 


37 


Summary 


SUBROUTINE  TRINV 

- The  subroutine  obtains  polar  coordinates  from  (two- 
dimensional)  cartesian  coordinates. 


Language 

Author 

Subroutine  statement 
Input  arguments 
Y 
X 

Output  arguments 
R 

TH 

Use  of  COMMON 
Source  deck 
Local  storage  used 


- ASA  Fortran  (Standard  Fortran  4). 

- R.H.  Gooding  (May  1968). 

- SUBROUTINE  TRINV  (Y,  X,  R,  TH) . 

Cartesian  y-coordinate  (arbitrary). 

Cartesian  x-coordinate  (arbitrary). 

Polar  r-coordinate . 

Polar  9-coordinate. 

- None. 

- 8 cards,  including  2 comment  cards  (ICL  code). 

- None . 


Subordinate  subprograms  - None. 

Explanation  - The  equations  r cos  9 = x , and  r sin  9 = y are 

solved  for  r and  9 , using  the  standard  ATAN2  function  to  give  0 between 
-7r  and  +tt  . If  x = y = 0 the  ATAN2  function  is  not  used  and  0 is  set  to 
zero . 


Remarks 

(1)  The  subroutine,  if  used  twice,  provides  a convenient  solution  of  the  three- 
dimensional  car tesian-to-po lar  transformation.  Thus  to  solve  the  equations. 


and 


r sin  9 cos  = x , 
r sin  0 sin  <P  - y , 
r cos  0 = z 


for  r,  0 and 
and 


the  following  two  statements  will  suffice: 

CALL  TRINV  (Y,  X,  RSINTH,  PHI)  , giving  r sin  0 and  <f> 

CALL  TRINV  (RSINTH,  Z,  R,  TH)  , giving  r and  0 . 


88 


Appendix  B 


Moreover,  this  solution  will  give  maximum  accuracy  for  0 and  <j>  , with  both 
angles  set  in  the  correct  quadrant. 

(2)  The  actual  input  and  output  arguments  must  be  distinct. 
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No. 

1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Author 

G.E.  Cook 
M.D.  Palmer 

R.H.  Merson 

P.R.  Peabody 
J.F.  Scott 
E.G.  Orozco 

P.R.  Peabody 
J.F.  Scott 
E.G.  Orozco 

G.E.  Cook 

R.H.  Gooding 
R.J.  Taylor 

Y.  Kozai 


Y.  Kozai 


A.  Cazenave 
F.  Forestier 

F.  Nouel 
J.L.  Pieplu 

D.G.  King-Hele 

G. E.  Cook 
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Fig. 4 Specimen  data  decks  for  POINT 
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